SystMaker.cxx
Go to the documentation of this file.
2 
6 
11 
12 #include "TCanvas.h"
13 #include "TLegend.h"
14 #include "TLegendEntry.h"
15 #include "TObjString.h"
16 #include "TKey.h"
17 
18 using std::cout;
19 using std::cerr;
20 using std::endl;
21 using std::string;
22 using std::vector;
23 using std::map;
24 using std::pair;
25 using std::make_pair;
26 using std::find;
27 using std::move;
28 using std::unique_ptr;
29 using std::make_unique;
30 using std::showpos;
31 using std::ostream;
32 using std::ostringstream;
33 using std::runtime_error;
34 
35 using osc::IOscCalc;
36 
37 using ana::covmx::Sample;
38 
39 namespace ana {
40 
41  // -------------------------------------------------------------------------------------
43  const HistAxis* axis, const HistAxis* numuAxis,
44  const Cut* FDCut, const Cut* NDCut, const Cut* numuCut,
45  const SystShifts* shiftData, const SystShifts* shiftMC,
46  const Var* weight,
47  Loaders* loaders, Loaders* loadersND) {
48 
49  IPrediction* pred = nullptr;
50 
51  if (predType == SystPredType::kExtrap) {
52  SterileGenerator genExtrap(*axis, *numuAxis,
53  *FDCut, *NDCut, *numuCut,
54  *shiftData, *weight);
55  if (loadersND)
56  pred = genExtrap.Generate(*loadersND, *loaders, *shiftMC).release();
57  else
58  pred = genExtrap.Generate(*loaders, *shiftMC).release();
59  }
60 
61  else if (predType == SystPredType::kFD) {
62  FDPredictionGenerator genFD(*axis, *FDCut, *shiftData, *weight);
63  pred = genFD.Generate(*loaders, *shiftMC).release();
64  }
65 
66  else if (predType == SystPredType::kND) {
67  NDPredictionGenerator genND(*axis, *NDCut, *shiftData, *weight);
68  pred = genND.Generate(*loaders, *shiftMC).release();
69  }
70 
71  return pred;
72 
73  } // function GetPrediction
74 
75  // --------------------------------------------------------------------------
76  TH1* GetSpectrum(IPrediction* pred, bool signal, const Sample sample,
77  IOscCalc* calc) {
78 
79  // If NC selection, signal is NC and background is CC. Nice and easy.
80  if (sample.selection == covmx::Selection::kNC) {
81  return pred->PredictComponent(calc, Flavors::kAll,
82  signal ? Current::kNC : Current::kCC, Sign::kBoth).ToTH1(sample.GetPOT());
83  }
84 
85  // Get beam polarity and CC signal channel for this sample
86  Sign::Sign_t sign = (sample.polarity==covmx::kFHC) ? Sign::kNu : Sign::kAntiNu;
87  Flavors::Flavors_t flavor;
88  if (sample.selection == covmx::kCCNumu)
89  flavor = Flavors::kNuMuToNuMu;
90  else if (sample.selection == covmx::kCCNue)
91  flavor = Flavors::kNuMuToNuE;
92  else assert(false and "Selection not recognised!");
93 
94  // Signal is the specified CC channel
95  if (signal) {
96  return pred->PredictComponent(calc, flavor, Current::kCC,
97  sign).ToTH1(sample.GetPOT());
98  }
99  // Background is everything else
100  else {
101  // Start with correct sign backgrounds
102  Spectrum ret = pred->PredictComponent(calc, Flavors::Flavors_t(~flavor),
103  Current::kCC, sign);
104  // CC wrong sign backgrounds
105  ret += pred->PredictComponent(calc, Flavors::kAll,
106  Current::kCC, Sign::Sign_t(~sign));
107  // NC everything
108  ret += pred->PredictComponent(calc, Flavors::kAll, Current::kNC,
109  Sign::kBoth);
110  return ret.ToTH1(sample.GetPOT());
111  }
112 
113  } // function GetSpectrum
114 
115  // --------------------------------------------------------------------------
116  string MakeLatexCommandName(const string& instring){
117 
118  TString mystring(instring.c_str());
119 
120  vector<TString> in = {"0","1","2","3","4","5",
121  "6","7","8","9","_","#"," ",".","-","(",")"};
122  vector<TString> out = {"zero","one","two","three","four","five",
123  "six","seven","eight","nine","","","","","XX","",""};
124  for(size_t i = 0; i < in.size(); ++i)
125  mystring.ReplaceAll(in[i],out[i]);
126 
127  return string(mystring.Data());
128  }
129 
130  // -------------------------------------------------------------------------------------
131  /// SystMakerShiftSigma constructor
133  fName(name),
134  fSigma(sigma) {
135  }
136 
137  // -------------------------------------------------------------------------------------
138  /// SystMakerShiftSigma destructor
140  if (fPrediction) delete fPrediction;
141  if (fShiftedSpectrum.first) delete fShiftedSpectrum.first;
142  if (fShiftedSpectrum.second) delete fShiftedSpectrum.second;
143  if (fShiftedRatio.first) delete fShiftedRatio.first;
144  if (fShiftedRatio.second) delete fShiftedRatio.second;
145  }
146 
147  // -------------------------------------------------------------------------------------
148  /// Print the sigmas for a given shift
150  os << " " << showpos << fSigma << ": "
151  << fShift.first << " (signal), "
152  << fShift.second << " (background)"
153  << endl;
154  }
155 
156  // -------------------------------------------------------------------------------------
157  void SystMakerShiftSigma::SaveTo(TDirectory* dir) const {
158 
159  TDirectory* tmp = gDirectory;
160 
161  dir->cd();
162  TObjString("SystMakerShiftSigma").Write("type");
163 
164  TObjString(fName.c_str()).Write("name");
165  TVectorD sigma(1);
166  sigma[0] = fSigma;
167  sigma.Write("sigma");
168 
169  fPrediction->SaveTo(dir, "pred");
170  TH1* shiftSig = fShiftedSpectrum.first;
171  TH1* shiftBkg = fShiftedSpectrum.second;
172  if (shiftSig)
173  shiftSig->Write("shiftSig");
174  if (shiftBkg)
175  shiftBkg->Write("shiftBkg");
176  TH1* shiftRatSig = fShiftedRatio.first;
177  TH1* shiftRatBkg = fShiftedRatio.second;
178  if (shiftRatSig)
179  shiftRatSig->Write("shiftRatSig");
180  if (shiftRatBkg)
181  shiftRatBkg->Write("shiftRatBkg");
182 
183  TVectorD v(2);
184  v[0] = fShift.first;
185  v[1] = fShift.second;
186  v.Write("shift");
187 
188  tmp->cd();
189 
190  }
191 
192  // -------------------------------------------------------------------------------------
194 
195  DontAddDirectory guard;
196 
197  TObjString* tag = (TObjString*)dir->Get("type");
198  assert(tag);
199  assert(tag->GetString() == "SystMakerShiftSigma");
200 
201  TObjString* name = (TObjString*)dir->Get("name");;
202  const TVectorD sigma = *(TVectorD*)dir->Get("sigma");
203  SystMakerShiftSigma* systMakerShiftSigma
204  = new SystMakerShiftSigma(name->GetString().Data(), sigma[0]);
205 
206  tag = (TObjString*)dir->Get("pred/type");
207  assert(tag);
208  if (tag->GetString() == "NDPredictionSterile")
209  systMakerShiftSigma->fPrediction = NDPredictionSterile::LoadFrom(dir, "pred").release();
210  else if (tag->GetString() == "FDPredictionSterile")
211  systMakerShiftSigma->fPrediction = FDPredictionSterile::LoadFrom(dir, "pred").release();
212  else if (tag->GetString() == "PredictionSterile")
213  systMakerShiftSigma->fPrediction = PredictionSterile::LoadFrom(dir, "pred").release();
214  else assert(false && "Prediction type not recognised!");
215 
216  TH1* shiftSig = (TH1*)dir->Get("shiftSig");
217  TH1* shiftBkg = (TH1*)dir->Get("shiftBkg");
218  if (shiftSig and shiftBkg)
219  systMakerShiftSigma->fShiftedSpectrum = make_pair(shiftSig, shiftBkg);
220  TH1* shiftRatSig = (TH1*)dir->Get("shiftRatSig");
221  TH1* shiftRatBkg = (TH1*)dir->Get("shiftRatBkg");
222  if (shiftRatSig and shiftRatBkg)
223  systMakerShiftSigma->fShiftedRatio = make_pair(shiftRatSig, shiftRatBkg);
224 
225  const TVectorD v = *(TVectorD*)dir->Get("shift");
226  systMakerShiftSigma->fShift = make_pair(v[0], v[1]);
227 
228  return systMakerShiftSigma;
229 
230  }
231 
232  // -------------------------------------------------------------------------------------
233  SystMakerShift::SystMakerShift(string name, string label) :
234  fName(name),
235  fLabel(label)
236  {}
237 
238  // -------------------------------------------------------------------------------------
239  /// Default destructor
241  {
242  for (auto sigma : fSigmas) {
243  if (sigma.second) delete sigma.second;
244  }
245  }
246 
247  // -------------------------------------------------------------------------------------
249  vector<int> sigmas;
250  for (const auto& sigma : fSigmas)
251  sigmas.push_back(sigma.first);
252  return sigmas;
253  }
254 
255  // -------------------------------------------------------------------------------------
256  void SystMakerShift::PrintShift(ostream& os) {
257  os << endl << " " << fLabel << " shifts:" << endl;
258  for (const auto& sigma : fSigmas)
259  sigma.second->PrintSigma(os);
260  return;
261  }
262 
263  // -------------------------------------------------------------------------------------
265  // Require -1, +1 shifts
266  bool good = true;
267  if (!fSigmas.count(-1) or !fSigmas.count(+1))
268  good = false;
269  return good;
270  }
271 
272  // -------------------------------------------------------------------------------------
273  TCanvas* SystMakerShift::DrawSig(TH1* nominal) {
274 
275  // Drawing tools ------------
276  TCanvas* canv = new TCanvas("canv", "", 800, 800);
277  TPad* padSpecSigTotal = new TPad("padSpecSigTotal", "", 0., 0.375, 1., 1.);
278  TPad* padRatioSigTotal = new TPad("padRatioSigTotal", "", 0., 0., 1., 0.375);
279  TLegend* leg = new TLegend(0.5, 0.5, 0.88, 0.88);
280  map<int, TLegendEntry*> legEntries;
281 
282  // Draw nominal ------------
283  TH1* hNom = (TH1*)nominal->Clone("hNom");
284  for (size_t bin = 1; bin <= size_t(hNom->GetNbinsX()); ++bin)
285  hNom->SetBinContent(bin, hNom->GetBinContent(bin)/hNom->GetBinWidth(bin));
286  hNom->GetYaxis()->SetTitle("Events / GeV");
287  hNom->SetLineStyle(kSolid);
288  hNom->SetLineColor(kBlack);
289  hNom->SetLineWidth(2);
290  padSpecSigTotal->cd();
291  hNom->Draw("hist");
292 
293  // ratio
294  TH1* nomRatio = (TH1*)nominal->Clone("nomRatio");
295  nomRatio->GetYaxis()->SetTitle("Shifted / Nominal");
296  for (Int_t bin = 0; bin <= nomRatio->GetNbinsX(); ++bin)
297  nomRatio->SetBinContent(bin, 1.);
298  padRatioSigTotal->cd();
299  nomRatio->Draw("hist");
300  legEntries[0] = new TLegendEntry(nominal, "Nominal", "l");
301  double ymaxRatio = 0;
302 
303  // Draw shifts ------------
304  int it = 0;
305  for (map<int, SystMakerShiftSigma*>::const_iterator sigma = fSigmas.begin(); sigma != fSigmas.end(); ++sigma, ++it) {
306  if (sigma->second->fSigma == 0)
307  continue;
308  // Total ------------
309  TH1* sigmaSig = (TH1*)sigma->second->fShiftedSpectrum.first->Clone();
310  for (size_t bin = 1; bin <= size_t(sigmaSig->GetNbinsX()); ++bin)
311  sigmaSig->SetBinContent(bin, sigmaSig->GetBinContent(bin)/sigmaSig->GetBinWidth(bin));
312  sigmaSig->GetYaxis()->SetTitle("Events / GeV");
313  sigmaSig->SetLineStyle(kDashed);
314  sigmaSig->SetLineColor(abs(sigma->second->fSigma)+1);
315  sigmaSig->SetLineWidth(2);
316  padSpecSigTotal->cd();
317  sigmaSig->DrawClone("hist same");
318  TH1* sigmaRatioSig = (TH1*)sigmaSig->Clone("sigmaRatioSig");
319  sigmaRatioSig->Divide(hNom);
320  sigmaRatioSig->SetBinContent(0, 1);
321  for (size_t i = 1; i <= (size_t)sigmaRatioSig->GetNbinsX(); ++i) {
322  double val = sigmaRatioSig->GetBinContent(i);
323  if (val == 0) {
324  sigmaRatioSig->SetBinContent(i, 1);
325  val = 1;
326  }
327  if (fabs(val-1) > ymaxRatio)
328  ymaxRatio = fabs(val-1);
329  }
330  padRatioSigTotal->cd();
331  sigmaRatioSig->DrawClone("hist same");
332  legEntries[abs(sigma->second->fSigma)]
333  = new TLegendEntry(sigmaSig->Clone(), Form("+/-%d#sigma", sigma->second->fSigma), "l");
334  }
335 
336  // Draw legend ------------
337  for (const auto& legEntry : legEntries)
338  leg->AddEntry(legEntry.second->GetObject(),
339  legEntry.second->GetLabel(),
340  legEntry.second->GetOption());
341  padSpecSigTotal->cd();
342  leg->Draw();
343 
344  nomRatio->SetMinimum(1-(1.1*ymaxRatio));
345  nomRatio->SetMaximum(1+(1.1*ymaxRatio));
346 
347  // Draw all together ------------
348  canv->cd();
349  canv->Clear();
350  padSpecSigTotal->Draw();
351  padRatioSigTotal->Draw();
352 
353  return canv;
354 
355  } // function SystMakerShift::DrawSig
356 
357  // -------------------------------------------------------------------------------------
358  TCanvas* SystMakerShift::DrawBkg(TH1* nominal) {
359 
360  // Drawing tools ------------
361  TCanvas* canv = new TCanvas("canv", "", 800, 800);
362  TPad* padSpecBkgTotal = new TPad("padSpecBkgTotal", "", 0., 0.375, 1., 1.);
363  TPad* padRatioBkgTotal = new TPad("padRatioBkgTotal", "", 0., 0., 1., 0.375);
364  TLegend* leg = new TLegend(0.5, 0.5, 0.88, 0.88);
365  map<int, TLegendEntry*> legEntries;
366 
367  // Draw nominal ------------
368  TH1* hNom = (TH1*)nominal->Clone("hNom");
369  for (size_t bin = 1; bin <= size_t(hNom->GetNbinsX()); ++bin)
370  hNom->SetBinContent(bin, hNom->GetBinContent(bin)/hNom->GetBinWidth(bin));
371  hNom->GetYaxis()->SetTitle("Events / GeV");
372  nominal->SetLineStyle(kSolid);
373  nominal->SetLineColor(kBlack);
374  nominal->SetLineWidth(2);
375  padSpecBkgTotal->cd();
376  nominal->Draw("hist");
377 
378  // ratio
379  TH1* nomRatio = (TH1*)nominal->Clone("nomRatio");
380  nomRatio->GetYaxis()->SetTitle("Shifted / Nominal");
381  for (Int_t bin = 0; bin <= nomRatio->GetNbinsX(); ++bin)
382  nomRatio->SetBinContent(bin, 1.);
383  padRatioBkgTotal->cd();
384  nomRatio->Draw("hist");
385  legEntries[0] = new TLegendEntry(nominal, "Nominal", "l");
386  double ymaxRatio = 0;
387 
388  // Draw shifts ------------
389  int it = 0;
390  for (map<int, SystMakerShiftSigma*>::const_iterator sigma = fSigmas.begin();
391  sigma != fSigmas.end(); ++sigma, ++it) {
392  if (sigma->second->fSigma == 0)
393  continue;
394  // Total ------------
395  TH1* sigmaBkg = (TH1*)sigma->second->fShiftedSpectrum.second->Clone();
396  for (size_t bin = 1; bin <= size_t(sigmaBkg->GetNbinsX()); ++bin)
397  sigmaBkg->SetBinContent(bin, sigmaBkg->GetBinContent(bin)/sigmaBkg->GetBinWidth(bin));
398  sigmaBkg->GetYaxis()->SetTitle("Events / GeV");
399  sigmaBkg->SetLineStyle(kDashed);
400  sigmaBkg->SetLineColor(abs(sigma->second->fSigma)+1);
401  sigmaBkg->SetLineWidth(2);
402  padSpecBkgTotal->cd();
403  sigmaBkg->DrawClone("hist same");
404  TH1* sigmaRatioBkg = (TH1*)sigmaBkg->Clone("sigmaRatioBkg");
405  sigmaRatioBkg->Divide(hNom);
406  sigmaRatioBkg->SetBinContent(0, 1);
407  for (size_t i = 1; i <= (size_t)sigmaRatioBkg->GetNbinsX(); ++i) {
408  double val = sigmaRatioBkg->GetBinContent(i);
409  if (val == 0) {
410  sigmaRatioBkg->SetBinContent(i, 1);
411  val = 1;
412  }
413  if (fabs(val-1) > ymaxRatio)
414  ymaxRatio = fabs(val-1);
415  }
416  padRatioBkgTotal->cd();
417  sigmaRatioBkg->DrawClone("hist same");
418  legEntries[abs(sigma->second->fSigma)]
419  = new TLegendEntry(sigmaBkg->Clone(), Form("+/-%d#sigma",
420  sigma->second->fSigma), "l");
421  }
422 
423  // Draw legend ------------
424  for (const auto& legEntry : legEntries)
425  leg->AddEntry(legEntry.second->GetObject(),
426  legEntry.second->GetLabel(),
427  legEntry.second->GetOption());
428  padSpecBkgTotal->cd();
429  leg->Draw();
430 
431  nomRatio->SetMinimum(1-(1.1*ymaxRatio));
432  nomRatio->SetMaximum(1+(1.1*ymaxRatio));
433 
434  // Draw all together ------------
435  canv->cd();
436  canv->Clear();
437  padSpecBkgTotal->Draw();
438  padRatioBkgTotal->Draw();
439 
440  return canv;
441 
442  } // function SystMakerShift::DrawBkg
443 
444  // -------------------------------------------------------------------------------------
445  void SystMakerShift::DrawShift(TDirectory* outDir, pair<TH1*, TH1*>& nominal) {
446 
447  // Out directories ------------
448  string dirSigShift = Form("Signal/Shifts/%s", fName.c_str());
449  string dirBkgShift = Form("Background/Shifts/%s", fName.c_str());
450  string dirSigShift_sigmas = dirSigShift+"/Sigmas/";
451  string dirBkgShiftSigmas = dirBkgShift+"/Sigmas/";
452  outDir->mkdir(dirSigShift.c_str());
453  outDir->mkdir(dirBkgShift.c_str());
454  outDir->mkdir(dirSigShift_sigmas.c_str());
455  outDir->mkdir(dirBkgShiftSigmas.c_str());
456 
457  // Drawing tools ------------
458  TCanvas* canv = new TCanvas("canv", "", 800, 800);
459  TPad* padSpecSig = new TPad("padSpecSig", "", 0., 0.375, 1., 1.);
460  TPad* padRatioSig = new TPad("padRatioSig", "", 0., 0., 1., 0.375);
461  TPad* padSpec_bg = new TPad("padSpec_bg", "", 0., 0.375, 1., 1.);
462  TPad* padRatio_bg = new TPad("padRatio_bg", "", 0., 0., 1., 0.375);
463  TPad* padSpecSigTotal = new TPad("padSpecSigTotal", "", 0., 0.375, 1., 1.);
464  TPad* padRatioSigTotal = new TPad("padRatioSigTotal", "", 0., 0., 1., 0.375);
465  TPad* padSpecBkgTotal = new TPad("padSpecBkgTotal", "", 0., 0.375, 1., 1.);
466  TPad* padRatioBkgTotal = new TPad("padRatioBkgTotal", "", 0., 0., 1., 0.375);
467  TLegend* leg = new TLegend(0.5, 0.5, 0.88, 0.88);
468  map<int, TLegendEntry*> legEntries;
469 
470  // Draw nominal ------------
471  // signal
472  TH1* nomSig = (TH1*)nominal.first->Clone("nomSig");
473  for (size_t bin = 1; bin <= size_t(nomSig->GetNbinsX()); ++bin)
474  nomSig->SetBinContent(bin, nomSig->GetBinContent(bin)/nomSig->GetBinWidth(bin));
475  nomSig->GetYaxis()->SetTitle("Events / GeV");
476  nomSig->SetLineStyle(kSolid);
477  nomSig->SetLineColor(kBlack);
478  nomSig->SetLineWidth(2);
479  padSpecSigTotal->cd();
480  nomSig->Draw("hist");
481  // background
482  TH1* nomBkg = (TH1*)nominal.second->Clone("nomBkg");
483  for (size_t bin = 1; bin <= size_t(nomBkg->GetNbinsX()); ++bin)
484  nomBkg->SetBinContent(bin, nomBkg->GetBinContent(bin)/nomBkg->GetBinWidth(bin));
485  nomBkg->GetYaxis()->SetTitle("Events / GeV");
486  nomBkg->SetLineStyle(kSolid);
487  nomBkg->SetLineColor(kBlack);
488  nomBkg->SetLineWidth(2);
489  padSpecBkgTotal->cd();
490  nomBkg->Draw("hist");
491  // ratio
492  TH1* nomRatio = (TH1*)nomSig->Clone("nomRatio");
493  nomRatio->GetYaxis()->SetTitle("Shifted / Nominal");
494  for (Int_t bin = 1; bin <= nomRatio->GetNbinsX(); ++bin)
495  nomRatio->SetBinContent(bin, 1.);
496  padRatioSigTotal->cd();
497  nomRatio->Draw("hist");
498  padRatioBkgTotal->cd();
499  nomRatio->Draw("hist");
500  legEntries[0] = new TLegendEntry(nomSig, "Nominal", "l");
501 
502  // Draw shifts ------------
503  int it = 0;
504  for (map<int, SystMakerShiftSigma*>::const_iterator sigma = fSigmas.begin(); sigma != fSigmas.end(); ++sigma, ++it) {
505  if (sigma->second->fSigma == 0)
506  continue;
507  // Total ------------
508  // signal
509  TH1* sigmaSig = (TH1*)sigma->second->fShiftedSpectrum.first->Clone();
510  for (size_t bin = 1; bin <= size_t(sigmaSig->GetNbinsX()); ++bin)
511  sigmaSig->SetBinContent(bin, sigmaSig->GetBinContent(bin)/sigmaSig->GetBinWidth(bin));
512  sigmaSig->GetYaxis()->SetTitle("Events / GeV");
513  sigmaSig->SetLineStyle(kDashed);
514  sigmaSig->SetLineColor(abs(sigma->second->fSigma)+1);
515  sigmaSig->SetLineWidth(2);
516  padSpecSigTotal->cd();
517  sigmaSig->DrawClone("hist same");
518  TH1* sigmaRatioSig = (TH1*)sigmaSig->Clone("sigmaRatioSig");
519  sigmaRatioSig->Divide(nomSig);
520  padRatioSigTotal->cd();
521  sigmaRatioSig->DrawClone("hist same");
522  // background
523  TH1* sigmaBkg = (TH1*)sigma->second->fShiftedSpectrum.second->Clone();
524  for (size_t bin = 1; bin <= size_t(sigmaBkg->GetNbinsX()); ++bin)
525  sigmaBkg->SetBinContent(bin, sigmaBkg->GetBinContent(bin)/sigmaBkg->GetBinWidth(bin));
526  sigmaBkg->GetYaxis()->SetTitle("Events / GeV");
527  sigmaBkg->SetLineStyle(kDashed);
528  sigmaBkg->SetLineColor(abs(sigma->second->fSigma)+1);
529  sigmaBkg->SetLineWidth(2);
530  padSpecBkgTotal->cd();
531  sigmaBkg->DrawClone("hist same");
532  TH1* sigmaRatioBkg = (TH1*)sigmaBkg->Clone("sigmaRatioBkg");
533  sigmaRatioBkg->Divide(nomBkg);
534  padRatioBkgTotal->cd();
535  sigmaRatioBkg->DrawClone("hist same");
536  legEntries[abs(sigma->second->fSigma)]
537  = new TLegendEntry(sigmaSig->Clone(), Form("+/-%d#sigma", sigma->second->fSigma), "l");
538  // Individual ------------
539  // signal
540  padSpecSig->cd();
541  padSpecSig->Clear();
542  nomSig->Draw("hist");
543  sigmaSig->SetLineColor(kRed);
544  sigmaSig->Draw("hist same");
545  padRatioSig->cd();
546  padRatioSig->Clear();
547  nomRatio->Draw("hist");
548  sigmaRatioSig->SetLineColor(kRed);
549  sigmaRatioSig->Draw("hist same");
550  canv->cd();
551  canv->Clear();
552  padSpecSig->DrawClone();
553  padRatioSig->DrawClone();
554  outDir->cd(dirSigShift_sigmas.c_str());
555  canv->Write(Form("%s%dSigma", fName.c_str(), sigma->second->fSigma));
556  // background
557  padSpec_bg->cd();
558  padSpec_bg->Clear();
559  nomBkg->Draw("hist");
560  sigmaBkg->SetLineColor(kRed);
561  sigmaBkg->Draw("hist same");
562  padRatio_bg->cd();
563  padRatio_bg->Clear();
564  nomRatio->Draw("hist");
565  sigmaRatioBkg->SetLineColor(kRed);
566  sigmaRatioBkg->Draw("hist same");
567  canv->cd();
568  canv->Clear();
569  padSpec_bg->DrawClone();
570  padRatio_bg->DrawClone();
571  outDir->cd(dirBkgShiftSigmas.c_str());
572  canv->Write(Form("%s%dSigma", fName.c_str(), sigma->second->fSigma));
573 
574  }
575 
576  // Draw legend ------------
577  for (const auto& legEntry : legEntries)
578  leg->AddEntry(legEntry.second->GetObject(),
579  legEntry.second->GetLabel(),
580  legEntry.second->GetOption());
581  padSpecSigTotal->cd();
582  leg->Draw();
583  padSpecBkgTotal->cd();
584  leg->Draw();
585 
586  // Draw all together ------------
587  // signal
588  canv->cd();
589  canv->Clear();
590  padSpecSigTotal->Draw();
591  padRatioSigTotal->Draw();
592  outDir->cd(dirSigShift.c_str());
593  canv->Write(Form("%sSigAllSigmas", fName.c_str()));
594  // background
595  canv->cd();
596  canv->Clear();
597  padSpecBkgTotal->Draw();
598  padRatioBkgTotal->Draw();
599  outDir->cd(dirBkgShift.c_str());
600  canv->Write(Form("%sBgAllSigmas", fName.c_str()));
601 
602  // clean up
603  // need to figure out what needs to be deleted
604  delete canv;
605 
606  } // function SystMakerShift::DrawShift
607 
608  // --------------------------------------------------------------------------
609  string SystMakerShift::WriteTable(pair<TH1*, TH1*> nominal, string name) {
610 
611  ostringstream table;
612 
613  table << "\\newcommand{\\" << name << "}\n";
614  table << "{\\centerline{\n";
615  table << "\\begin{tabular}{ l | r r r | r r }\n";
616  table << " & Nom.&(+)&(-)&\\%(+)&\\%(-) \\\\ \\hline\n";
617 
618  // We write +1 and -1 sigmas to a table. First check they both exist!
619 
620  vector<int> sigmas = this->GetSigmas();
621  if (find(sigmas.begin(), sigmas.end(), 1) == sigmas.end())
622  assert(false and "+1 shift required to call WriteTable()!");
623  bool useDown = find(sigmas.begin(), sigmas.end(), -1) != sigmas.end();
624 
625  // Get nominal event counts for signal and background
626  double sigNom = 0, bkgNom = 0;
627  size_t nbins = nominal.first->GetNbinsX();
628  for (size_t it = 1; it <= nbins; ++it) {
629  sigNom += nominal.first->GetBinContent(it);
630  bkgNom += nominal.second->GetBinContent(it);
631  }
632 
633  // Now calculate the up and down totals
634  double sigUp = sigNom, bkgUp = bkgNom, sigDown = 0, bkgDown = 0;
635  if (useDown) {
636  sigDown = sigNom;
637  bkgDown = bkgNom;
638  }
639 
640  for (size_t it = 1; it <= nbins; ++it) {
641  sigUp += this->GetShiftedSpectrum(1).first->GetBinContent(it)
642  - nominal.first->GetBinContent(it);
643  sigDown += this->GetShiftedSpectrum(-1).first->GetBinContent(it)
644  - nominal.first->GetBinContent(it);
645  bkgUp += this->GetShiftedSpectrum(1).second->GetBinContent(it)
646  - nominal.second->GetBinContent(it);
647  if (useDown) {
648  sigDown += this->GetShiftedSpectrum(-1).first->GetBinContent(it)
649  - nominal.first->GetBinContent(it);
650  bkgDown += this->GetShiftedSpectrum(-1).second->GetBinContent(it)
651  - nominal.second->GetBinContent(it);
652  }
653  }
654 
655  // Write them to the table and return it
656  table << "Signal & " << sigNom << " & " << sigUp << " & ";
657  if (useDown) table << sigDown;
658  else table << "N/A";
659  table << " & " << 100*((sigUp/sigNom)-1) << " & ";
660  if (useDown) table << 100*((sigDown/sigNom)-1);
661  else table << "N/A";
662  table << " \\\\\n";
663  table << "Background & " << bkgNom << " & " << bkgUp << " & ";
664  if (useDown) table << bkgDown;
665  else table << "N/A";
666  table << " & " << 100*((bkgUp/bkgNom)-1) << " & ";
667  if (useDown) table << 100*((bkgDown/bkgNom)-1);
668  else table << "N/A";
669  table << " \\\\\n\\end{tabular}\n}}\n";
670 
671  return table.str();
672 
673  } // function SystMakerShift::WriteTable
674 
675  // -------------------------------------------------------------------------------------
676  pair<TH1*, TH1*> SystMakerShift::GetShiftedSpectrum(int sigma) {
677  if (!fSigmas.count(sigma)) {
678  cout << "Error: attempting to retrieve " << sigma << " sigma shift from "
679  << fName << "; does not exist. Aborting..." << endl;
680  abort();
681  }
682  return fSigmas[sigma]->fShiftedSpectrum;
683  }
684 
685  // -------------------------------------------------------------------------------------
686  void SystMakerShift::ProcessShift(pair<TH1*, TH1*>& nominal_spectrum,
687  const Sample sample, IOscCalc* calc) {
688  for (auto& sigma : fSigmas) {
689  TH1* hShiftSig = GetSpectrum(sigma.second->fPrediction, true, sample, calc);
690  TH1* hShiftBg = GetSpectrum(sigma.second->fPrediction, false, sample, calc);
691  TH1* hShiftRatioSig = (TH1*)hShiftSig->Clone("shiftRatSig");
692  hShiftRatioSig->Divide(nominal_spectrum.first);
693  TH1* hShiftRatioBg = (TH1*)hShiftBg->Clone("shiftRatBkg");
694  hShiftRatioBg->Divide(nominal_spectrum.second);
695  double shiftSig = abs(hShiftSig->Integral() - nominal_spectrum.first->Integral()) /
696  nominal_spectrum.first->Integral();
697  double shiftBg = abs(hShiftBg->Integral() - nominal_spectrum.second->Integral()) /
698  nominal_spectrum.second->Integral();
699  sigma.second->fShiftedSpectrum = make_pair(hShiftSig, hShiftBg);
700  sigma.second->fShiftedRatio = make_pair(hShiftRatioSig, hShiftRatioBg);
701  sigma.second->fShift = make_pair(shiftSig, shiftBg);
702  }
703  return;
704  }
705 
706  // -------------------------------------------------------------------------------------
707  /// Base class implementation of MakePredictions, which just throws an error
709  const HistAxis* axis, const HistAxis* numuAxis,
710  const Cut* FDCut, const Cut* NDCut, const Cut* numuCut,
711  const Var* weight, Loaders* loaders) {
712 
713  // This function has to be defined in the base class for <reasons> but the
714  // user should never actually call it so it just throws an error
715 
716  assert(false and "Base class version of MakePredictions should never be called!");
717 
718  } // function SystMakerShift::MakePredictions
719 
720  // -------------------------------------------------------------------------------------
721  /// SaveTo implementation for SystMakerShift
722  void SystMakerShift::SaveTo(TDirectory* dir) const {
723 
724  TDirectory* tmp = gDirectory;
725 
726  dir->cd();
727  TObjString("SystMakerShift").Write("type");
728 
729  TObjString(fName.c_str()).Write("name");
730  TObjString(fLabel.c_str()).Write("label");
731 
732  const char* sigmaDir = "Sigmas";
733  dir->mkdir(sigmaDir);
734  for (const auto& sigma : fSigmas) {
735  const char* sigmaSubdir = Form("%s/Sigma%d", sigmaDir, sigma.second->fSigma);
736  dir->mkdir(sigmaSubdir);
737  sigma.second->SaveTo(dir->GetDirectory(sigmaSubdir));
738  }
739 
740  tmp->cd();
741 
742  } // function SystMakerShift::SaveTo
743 
744  // -------------------------------------------------------------------------------------
745  /// LoadFrom implementation for SystMakerShift
746  unique_ptr<SystMakerShift> SystMakerShift::LoadFrom(TDirectory* dir) {
747 
748  DontAddDirectory guard;
749 
750  TObjString* tag = (TObjString*)dir->Get("type");
751  assert(tag);
752  assert(tag->GetString() == "SystMakerShift");
753 
754  TObjString* name = (TObjString*)dir->Get("name");;
755  TObjString* label = (TObjString*)dir->Get("label");;
756  unique_ptr<SystMakerShift> systMakerShift
757  = make_unique<SystMakerShift>(name->GetString().Data(), label->GetString().Data());
758 
759  map<int, SystMakerShiftSigma*> sigmas;
760  TIter next(dir->GetDirectory("Sigmas")->GetListOfKeys());
761  TKey* key;
762  while ((key = (TKey*)next())) {
763  SystMakerShiftSigma* sigma
764  = SystMakerShiftSigma::LoadFrom(dir->GetDirectory(Form("Sigmas/%s", key->GetName())));
765  sigmas[sigma->fSigma] = sigma;
766  }
767  systMakerShift->fSigmas = sigmas;
768 
769  return systMakerShift;
770 
771  } // function SystMakerShift::LoadFrom
772 
773  // -------------------------------------------------------------------------------------
774  /// SystMakerLoaderShift constructor
776  SystMakerShift(name, label) {
777  }
778 
779  // -------------------------------------------------------------------------------------
780  /// Add sigma for loader shift with the same ND and FD loader
782  int sigma) {
783 
784  this->AddSigma(loaders, loaders, sigma);
785 
786  } // function SystMakerLoaderShift::AddSigma
787 
788  // -------------------------------------------------------------------------------------
789  /// Add sigma for loader shift with separate ND and FD loaders
791  Loaders* loaders_fd,
792  int sigma) {
793 
794  SystMakerShiftSigma* shiftSigma = new SystMakerShiftSigma(fName, sigma);
795  fSigmas[sigma] = shiftSigma;
796  fLoaders[sigma] = make_pair(loadersND, loaders_fd);
797 
798  } // function SystMakerLoaderShift::AddSigma
799 
800  // -------------------------------------------------------------------------------------
801  /// Add on/off for loader shift
803 
804  this->AddOnOff(loaders, loaders);
805 
806  } // function SystMakerLoaderShift::AddOnOff
807 
808  // -------------------------------------------------------------------------------------
809  /// Add on/off for loader shift with separate ND and FD loaders
811  Loaders* loaders_fd) {
812 
813  for (const int sigma : {-1,1}) {
814  SystMakerShiftSigma* shiftSigma = new SystMakerShiftSigma(fName, sigma);
815  fSigmas[sigma] = shiftSigma;
816  fLoaders[sigma] = pair<Loaders*, Loaders*>(loadersND, loaders_fd);
817  }
818 
819  // std::cout << "Calling CheckLoaders in AddOnOff" << std::endl;
820  // CheckLoaders();
821 
822  } // function SystMakerLoaderShift::AddOnOff
823 
824  // -------------------------------------------------------------------------------------
825  /// Make shifted prediction for systematic shift derived from alternate MC loader
827  const HistAxis* axis, const HistAxis* numuAxis,
828  const Cut* FDCut, const Cut* NDCut, const Cut* numuCut,
829  const Var* weight, Loaders* loaders) {
830 
831  // CheckLoaders();
832 
833  for (auto& sigma : fSigmas) {
834 
835  sigma.second->fPrediction = GetPrediction(predType,
836  axis, numuAxis,
837  FDCut, NDCut, numuCut,
838  &kNoShift, &kNoShift,
839  weight,
840  fLoaders[sigma.second->fSigma].second,
841  fLoaders[sigma.second->fSigma].first);
842  }
843 
844  } // function SystMakerLoaderShift::MakePredictions
845 
846  // -------------------------------------------------------------------------------------
847  /// SystMakerWeightShift constructor
849  SystMakerShift(name, label) {
850  }
851 
852  // -------------------------------------------------------------------------------------
853  /// SystMakerWeightShift constructor
855  for (auto& weight : fWeights)
856  delete weight.second;
857  }
858 
859  // -------------------------------------------------------------------------------------
860  /// Add sigma for weight shift
861  void SystMakerWeightShift::AddSigma(const Var* weight, int sigma) {
862 
863  SystMakerShiftSigma* shiftSigma = new SystMakerShiftSigma(fName, sigma);
864  fSigmas[sigma] = shiftSigma;
865  fWeights[sigma] = weight;
866  } // function SystMakerWeightShift::AddSigma
867 
868  // -------------------------------------------------------------------------------------
869  /// Add on/off for weight shift
871 
872  for (const auto& sigma : {-1,1}) {
873  SystMakerShiftSigma* shiftSigma = new SystMakerShiftSigma(fName, sigma);
874  fSigmas[sigma] = shiftSigma;
875  fWeights[sigma] = weight;
876  }
877 
878  } // function SystMakerWeightShift::AddOnOff
879 
880  // -------------------------------------------------------------------------------------
881  /// Make shifted prediction for systematic shift derived from CAFAna weight
883  const HistAxis* axis, const HistAxis* numuAxis,
884  const Cut* FDCut, const Cut* NDCut, const Cut* numuCut,
885  const Var* weight, Loaders* loaders) {
886 
887  for (auto& sigma : fSigmas)
888  sigma.second->fPrediction = GetPrediction(predType, axis, numuAxis, FDCut,
889  NDCut, numuCut, &kNoShift, &kNoShift, fWeights[sigma.second->fSigma], loaders);
890 
891  } // function SystMakerWeightShift::MakePrediction
892 
893  // -------------------------------------------------------------------------------------
894  /// SystMakerSystShift constructor
896  SystMakerShift(name, label) {
897  }
898 
899  // -------------------------------------------------------------------------------------
900  /// Add sigma for syst shift
901  void SystMakerSystShift::AddSigma(SystShifts* syst, int sigma) {
902  SystMakerShiftSigma* shiftSigma = new SystMakerShiftSigma(fName, sigma);
903  fSigmas[sigma] = shiftSigma;
904  fSysts[sigma] = syst;
905  return;
906  }
907 
908  // -------------------------------------------------------------------------------------
909  /// Make shifted prediction for systematic shift derived from SystShift
911  const HistAxis* axis, const HistAxis* numuAxis,
912  const Cut* FDCut, const Cut* NDCut, const Cut* numuCut,
913  const Var* weight, Loaders* loaders) {
914 
915  for (auto& sigma : fSigmas) {
916  sigma.second->fPrediction = GetPrediction(predType, axis, numuAxis, FDCut,
917  NDCut, numuCut, &kNoShift, fSysts[sigma.second->fSigma], weight, loaders);
918  }
919  } // function SystMakerSystShift::MakePrediction
920 
921  // -------------------------------------------------------------------------------------
922  /// SystMaker constructor
923  SystMaker::SystMaker(string name, string label, string sampleName) :
924  fName(name),
925  fLabel(label),
926  fSampleName(sampleName),
927  fAxis(nullptr),
928  fNuMuAxis(nullptr),
929  fNDCut(nullptr),
930  fFDCut(nullptr),
931  fNuMuCut(nullptr),
932  fLoaders(nullptr),
933  fWeight(nullptr),
934  fNominalPred(nullptr) {
935 
936  } // SystMaker constructor
937 
938  // -------------------------------------------------------------------------------------
939  /// SystMaker destructor
941 
942  for (auto shift : fShifts) {
943  if (shift.second) delete shift.second;
944  }
945  if (fNominalPred) delete fNominalPred;
946 
947  } // SystMaker destructor
948 
949  // -------------------------------------------------------------------------------------
950  /// Get the shift for a given name
951  SystMakerShift* SystMaker::GetShift(string shift_name) {
952  if (!fShifts.count(shift_name)) {
953  cout << "Error: requested shift " << shift_name << " does not exist. Aborting..." << endl;
954  abort();
955  }
956  return fShifts[shift_name];
957 
958  } // function SystMaker::GetShift
959 
960  // -------------------------------------------------------------------------------------
961  /// Get the total shifts for a given sigma (sig+bkg)
962  pair<double, double> SystMaker::GetTotalShift(int sigma) {
963  if (!fTotalShifts.count(sigma)) {
964  cout << "Error: sigma " << sigma << " does not exist. Aborting..." << endl;
965  abort();
966  }
967  return fTotalShifts[sigma];
968 
969  } // function SystMakerShift::GetTotalShift
970 
971  // -------------------------------------------------------------------------------------
972  /// Get the total shift for a given sigma and channel selection (sig/bkg)
973  double SystMaker::GetTotalShift(int sigma, bool signal) {
974 
975  pair<double, double> totalShift = this->GetTotalShift(sigma);
976  return signal ? totalShift.first : totalShift.second;
977 
978  } // function SystMakerShift::GetTotalShift
979 
980  // -------------------------------------------------------------------------------------
981  /// Get the sigmas for a given systematic shift
982  vector<int> SystMaker::GetSigmas() {
983 
984  map<int, vector<string> > shifts_sigmas;
985  for (const auto& shift : fShifts) {
986  vector<int> shiftSigmas = shift.second->GetSigmas();
987  for (const auto& shiftSigma : shiftSigmas)
988  shifts_sigmas[shiftSigma].push_back(shift.second->GetName());
989  }
990  vector<int> sigmas;
991  for (const auto& shiftSigmas : shifts_sigmas)
992  if (shiftSigmas.second.size() == fShifts.size())
993  sigmas.push_back(shiftSigmas.first);
994 
995  return sigmas;
996 
997  } // function SystMaker::GetSigmas
998 
999  // -------------------------------------------------------------------------------------
1000  /// Print the total shifts for each shift owned by this syst maker
1001  void SystMaker::PrintSyst(ostream& os) {
1002 
1003  os << endl << endl << showpos
1004  << "Systematic " << fLabel << " has total shift:" << endl
1005  << " -1: " << fTotalShifts[-1].first << " (signal), " << fTotalShifts[-1].second << " (background)" << endl
1006  << " +1: " << fTotalShifts[+1].first << " (signal), " << fTotalShifts[+1].second << " (background)" << endl;
1007  os << endl << "Individual shifts:" << endl;
1008  for (const auto& shift : fShifts)
1009  shift.second->PrintShift(os);
1010 
1011  } // function SystMaker::PrintSyst
1012 
1013  // -------------------------------------------------------------------------------------
1014  /// Add SystMakerShift to this syst maker
1016 
1017  fShifts[shift->GetName()] = shift;
1018 
1019  } // function SystMaker::AddShift
1020 
1021  // -------------------------------------------------------------------------------------
1022  /// Make plots of all the shifted spectra and ratios for this systmaker
1023  TCanvas* SystMaker::DrawSig() {
1024 
1025  // Nominal
1026  TH1* nomSig = (TH1*)fSpectra[0].first->Clone();
1027  for (size_t bin = 1; bin <= size_t(nomSig->GetNbinsX()); ++bin)
1028  nomSig->SetBinContent(bin, nomSig->GetBinContent(bin)/nomSig->GetBinWidth(bin));
1029  nomSig->GetYaxis()->SetTitle("Events / GeV");
1030  nomSig->SetLineColor(kBlack);
1031  nomSig->SetLineStyle(kSolid);
1032  nomSig->SetLineWidth(2);
1033  TH1* nomRatio = (TH1*)nomSig->Clone();
1034  nomRatio->GetYaxis()->SetTitle("Shifted / Nominal");
1035  for (size_t bin = 0; bin <= (size_t)nomRatio->GetNbinsX(); ++bin)
1036  nomRatio->SetBinContent(bin, 1.);
1037 
1038  // Draw total shifts -- all sigmas
1039  TCanvas* canv = new TCanvas("canvSig", "", 800, 800);
1040  TPad* padSpec = new TPad("padSpecSig", "", 0., 0.375, 1., 1. );
1041  TPad* padRatio = new TPad("padRatioSig", "", 0., 0., 1., 0.375);
1042  map<int, TLegendEntry*> legEntrySigma;
1043  padSpec->cd();
1044  nomSig->Draw("hist");
1045  padRatio->cd();
1046  nomRatio->Draw("hist");
1047  legEntrySigma[0] = new TLegendEntry(nomSig, "Nominal", "l");
1048  double ymaxRatio = 0;
1049 
1050  // Loop over each sigma
1051  for (const auto& sigma : this->GetSigmas()) {
1052  if (sigma == 0)
1053  continue;
1054 
1055  // Draw the spectrum
1056  TH1* sigmaSig = (TH1*)fSpectra[sigma].first->Clone();
1057  for (size_t bin = 1; bin <= size_t(sigmaSig->GetNbinsX()); ++bin)
1058  sigmaSig->SetBinContent(bin, sigmaSig->GetBinContent(bin)/sigmaSig->GetBinWidth(bin));
1059  sigmaSig->GetYaxis()->SetTitle("Events / GeV");
1060  double hmax = nomSig->GetMaximum() > sigmaSig->GetMaximum() ?
1061  nomSig->GetMaximum() : sigmaSig->GetMaximum();
1062  nomSig->SetMaximum(1.1*hmax);
1063  sigmaSig->SetLineColor(abs(sigma)+1);
1064  sigmaSig->SetLineStyle(kDashed);
1065  sigmaSig->SetLineWidth(2);
1066  padSpec->cd();
1067  sigmaSig->Draw("hist same");
1068  if (!legEntrySigma.count(abs(sigma)))
1069  legEntrySigma[abs(sigma)] = new TLegendEntry(sigmaSig, Form("+/-%d #sigma", abs(sigma)), "l");
1070 
1071  // Draw the ratio
1072  TH1* ratioSig = (TH1*)sigmaSig->Clone();
1073  ratioSig->Divide(nomSig);
1074  ratioSig->SetBinContent(0, 1);
1075  for (size_t i = 1; i <= (size_t)ratioSig->GetNbinsX(); ++i) {
1076  double val = ratioSig->GetBinContent(i);
1077  if (val == 0) {
1078  ratioSig->SetBinContent(i, 1);
1079  val = 1;
1080  }
1081  if (fabs(1-val) > ymaxRatio)
1082  ymaxRatio = fabs(1-val);
1083  }
1084  padRatio->cd();
1085  ratioSig->Draw("hist same");
1086  }
1087  TLegend* legSigma = new TLegend(0.5, 0.5, 0.88, 0.88);
1088  for (const auto& legEntry : legEntrySigma)
1089  legSigma->AddEntry(legEntry.second->GetObject(),
1090  legEntry.second->GetLabel(),
1091  legEntry.second->GetOption());
1092  padSpec->cd();
1093  legSigma->Draw("same");
1094 
1095  nomRatio->SetMinimum(1-(1.1*ymaxRatio));
1096  nomRatio->SetMaximum(1+(1.1*ymaxRatio));
1097 
1098  canv->cd();
1099  canv->Clear();
1100  padSpec->Draw();
1101  padRatio->Draw();
1102 
1103  return canv;
1104 
1105  } // function SystMaker::DrawSig
1106 
1107  // -------------------------------------------------------------------------------------
1108  /// Make plots of all the shifted spectra and ratios for this systmaker
1109  TCanvas* SystMaker::DrawBkg() {
1110 
1111  // Nominal
1112  TH1* nomBkg = (TH1*)fSpectra[0].second->Clone();
1113  for (size_t bin = 1; bin <= size_t(nomBkg->GetNbinsX()); ++bin)
1114  nomBkg->SetBinContent(bin, nomBkg->GetBinContent(bin)/nomBkg->GetBinWidth(bin));
1115  nomBkg->GetYaxis()->SetTitle("Events / GeV");
1116  nomBkg->SetLineColor(kBlack);
1117  nomBkg->SetLineStyle(kSolid);
1118  nomBkg->SetLineWidth(2);
1119  TH1* nomRatio = (TH1*)nomBkg->Clone();
1120  nomRatio->GetYaxis()->SetTitle("Shifted / Nominal");
1121  for (size_t bin = 0; bin <= (size_t)nomRatio->GetNbinsX(); ++bin)
1122  nomRatio->SetBinContent(bin, 1.);
1123 
1124  // Draw total shifts -- all sigmas
1125  TCanvas* canv = new TCanvas("canvBkg", "", 800, 800);
1126  TPad* padSpec = new TPad("padSpecBkg", "", 0., 0.375, 1., 1. );
1127  TPad* padRatio = new TPad("padRatio_bkg", "", 0., 0., 1., 0.375);
1128  map<int, TLegendEntry*> legEntrySigma;
1129  padSpec->cd();
1130  nomBkg->Draw("hist");
1131  padRatio->cd();
1132  nomRatio->Draw("hist");
1133  legEntrySigma[0] = new TLegendEntry(nomBkg, "Nominal", "l");
1134  double ymaxRatio = 0;
1135 
1136  // Loop over each sigma
1137  for (const auto& sigma : this->GetSigmas()) {
1138  if (sigma == 0)
1139  continue;
1140 
1141  // Draw the spectrum
1142  TH1* sigmaBkg = (TH1*)fSpectra[sigma].second->Clone();
1143  for (size_t bin = 1; bin <= size_t(sigmaBkg->GetNbinsX()); ++bin)
1144  sigmaBkg->SetBinContent(bin, sigmaBkg->GetBinContent(bin)/sigmaBkg->GetBinWidth(bin));
1145  sigmaBkg->GetYaxis()->SetTitle("Events / GeV");
1146  double hmax = nomBkg->GetMaximum() > sigmaBkg->GetMaximum() ?
1147  nomBkg->GetMaximum() : sigmaBkg->GetMaximum();
1148  nomBkg->SetMaximum(1.1*hmax);
1149  sigmaBkg->SetLineColor(abs(sigma)+1);
1150  sigmaBkg->SetLineStyle(kDashed);
1151  sigmaBkg->SetLineWidth(2);
1152  padSpec->cd();
1153  sigmaBkg->Draw("hist same");
1154  if (!legEntrySigma.count(abs(sigma)))
1155  legEntrySigma[abs(sigma)] = new TLegendEntry(sigmaBkg, Form("+/-%d #sigma", abs(sigma)), "l");
1156 
1157  // Draw the ratio
1158  TH1* ratioBkg = (TH1*)sigmaBkg->Clone();
1159  ratioBkg->Divide(nomBkg);
1160  ratioBkg->SetBinContent(0, 1);
1161  for (size_t i = 1; i <= (size_t)ratioBkg->GetNbinsX(); ++i) {
1162  double val = ratioBkg->GetBinContent(i);
1163  if (val == 0) {
1164  ratioBkg->SetBinContent(i, 1);
1165  val = 1;
1166  }
1167  if (fabs(1-val) > ymaxRatio)
1168  ymaxRatio = fabs(1-val);
1169  }
1170  padRatio->cd();
1171  ratioBkg->Draw("hist same");
1172  }
1173  TLegend* legSigma = new TLegend(0.5, 0.5, 0.88, 0.88);
1174  for (const auto& legEntry : legEntrySigma)
1175  legSigma->AddEntry(legEntry.second->GetObject(),
1176  legEntry.second->GetLabel(),
1177  legEntry.second->GetOption());
1178  padSpec->cd();
1179  legSigma->Draw("same");
1180 
1181  nomRatio->SetMinimum(1-(1.1*ymaxRatio));
1182  nomRatio->SetMaximum(1+(1.1*ymaxRatio));
1183 
1184  canv->cd();
1185  canv->Clear();
1186  padSpec->Draw();
1187  padRatio->Draw();
1188 
1189  return canv;
1190 
1191  } // function SystMaker::DrawBkg
1192 
1193  // -------------------------------------------------------------------------------------
1194  /// Make plots of all the shifted spectra and ratios for this systmaker
1195  void SystMaker::DrawSyst(TDirectory* outDir) {
1196 
1197  string dirTop = fName;
1198  string dirSignal = dirTop+"/Signal/";
1199  string dirBackground = dirTop+"/Background/";
1200  string dirSigShifts = dirSignal+"/Shifts/";
1201  string dirBkgShifts = dirBackground+"/Shifts/";
1202  outDir->mkdir(dirTop.c_str());
1203  outDir->mkdir(dirSignal.c_str());
1204  outDir->mkdir(dirBackground.c_str());
1205  outDir->mkdir(dirSigShifts.c_str());
1206  outDir->mkdir(dirBkgShifts.c_str());
1207 
1208  TCanvas* canvSig = DrawSig();
1209  TCanvas* canvBkg = DrawBkg();
1210  outDir->cd(dirSignal.c_str());
1211  canvSig->Write("TotalShiftSig");
1212  outDir->cd(dirBackground.c_str());
1213  canvBkg->Write("TotalShiftBg");
1214 
1215  // Draw total shifts -- +1/-1 sigma
1216 
1217  // Draw individual shifts
1218  for (const auto& shift : fShifts)
1219  shift.second->DrawShift(outDir->GetDirectory(dirTop.c_str()), fSpectra[0]);
1220 
1221  // clean up
1222  delete canvSig;
1223  delete canvBkg;
1224 
1225  } // function SystMaker::DrawSyst
1226 
1227  // --------------------------------------------------------------------------
1228  /// Write a table with up & down shifts for this systematic
1229  string SystMaker::WriteTable(string name) {
1230 
1231  ostringstream table;
1232 
1233  table << "\\newcommand{\\" << name << "}\n";
1234  table << "{\\centerline{\n";
1235  table << "\\begin{tabular}{ l | r r r | r r }\n";
1236  table << " & Nom.&(+)&(-)&\\%(+)&\\%(-) \\\\ \\hline\n";
1237 
1238  // We write +1 and -1 sigmas to a table. First check they both exist!
1239  vector<int> sigmas = this->GetSigmas();
1240  if (find(sigmas.begin(), sigmas.end(), 1) == sigmas.end())
1241  assert(false and "+1 sigma shift required to call WriteTable()!");
1242  bool useDown = find(sigmas.begin(), sigmas.end(), -1) != sigmas.end();
1243 
1244  // Get nominal event counts for signal and background
1245  double sigNom = 0, bkgNom = 0;
1246  size_t nbins = fSpectra[0].first->GetNbinsX();
1247  for (size_t it = 1; it <= nbins; ++it) {
1248  sigNom += fSpectra[0].first->GetBinContent(it);
1249  bkgNom += fSpectra[0].second->GetBinContent(it);
1250  }
1251 
1252  // Now calculate the up and down totals
1253  double sigUp = sigNom, bkgUp = bkgNom, sigDown = 0, bkgDown = 0;
1254  if (useDown) {
1255  sigDown = sigNom;
1256  bkgDown = bkgNom;
1257  }
1258 
1259  for (auto& shift : fShifts) {
1260  for (size_t it = 1; it <= nbins; ++it) {
1261  sigUp += shift.second->GetShiftedSpectrum(1).first->GetBinContent(it)
1262  - fSpectra[0].first->GetBinContent(it);
1263  bkgUp += shift.second->GetShiftedSpectrum(1).second->GetBinContent(it)
1264  - fSpectra[0].second->GetBinContent(it);
1265  if (useDown) {
1266  sigDown += shift.second->GetShiftedSpectrum(-1).first->GetBinContent(it)
1267  - fSpectra[0].first->GetBinContent(it);
1268  bkgDown += shift.second->GetShiftedSpectrum(-1).second->GetBinContent(it)
1269  - fSpectra[0].second->GetBinContent(it);
1270  }
1271  }
1272  }
1273 
1274  // Write them to the table and return it
1275  table << "Signal & " << sigNom << " & " << sigUp << " & ";
1276  if (useDown) table << sigDown;
1277  else table << "N/A";
1278  table << " & " << 100*((sigUp/sigNom)-1) << " & ";
1279  if (useDown) table << 100*((sigDown/sigNom)-1);
1280  else table << "N/A";
1281  table << " \\\\\n";
1282  table << "Background & " << bkgNom << " & " << bkgUp << " & ";
1283  if (useDown) table << bkgDown;
1284  else table << "N/A";
1285  table << " & " << 100*((bkgUp/bkgNom)-1) << " & ";
1286  if (useDown) table << 100*((bkgDown/bkgNom)-1);
1287  else table << "N/A";
1288  table << " \\\\\n\\end{tabular}\n}}\n";
1289 
1290  return table.str();
1291 
1292  } // function SystMaker::WriteTable
1293 
1294  // -------------------------------------------------------------------------------------
1295  /// Make predictions for all the SystMakerShift objects owned by this syst maker
1297 
1298  fNominalPred = GetPrediction(predType,
1299  fAxis, fNuMuAxis,
1301  &kNoShift, &kNoShift,
1302  fWeight, fLoaders);
1303 
1304  for (const auto& shift : fShifts) {
1305  shift.second->MakePredictions(predType,
1306  fAxis, fNuMuAxis,
1308  fWeight, fLoaders);
1309  }
1310  } // function SystMaker::MakePredictions
1311 
1312  // -------------------------------------------------------------------------------------
1313  void SystMaker::ProcessSyst(const Sample sample, IOscCalc* calc) {
1314 
1315  // Get nominal spectrum
1316  fSpectra[0] = make_pair(GetSpectrum(fNominalPred, true, sample, calc),
1317  GetSpectrum(fNominalPred, false, sample, calc));
1318 
1319  // Process all shifts
1320  for (const auto& shift : fShifts) {
1321  shift.second->ProcessShift(fSpectra[0], sample, calc);
1322  SystMakerShift* systShift = static_cast<SystMakerShift*>(shift.second);
1323  fShifts[shift.first] = systShift;
1324  }
1325 
1326  // Get total +1/-1 sigma systematic shift -- required
1327  // Also get total sigma shift for all sigmas in this syst
1328  vector<int> sigmas = this->GetSigmas();
1329  Int_t nbins = fSpectra[0].first->GetNbinsX();
1330  for (const auto& sigma : sigmas) {
1331  vector<double> evtDiffSig(nbins, 0.), evtDiffBg(nbins, 0.);
1332  TH1* shiftSig = (TH1*)fSpectra[0].first->Clone("shiftSig");
1333  TH1* shiftBkg = (TH1*)fSpectra[0].second->Clone("shiftBkg");
1334  double totShiftSig = 0., totShiftBkg = 0.;
1335  for (Int_t bin = 1; bin <= nbins; ++bin) {
1336  double valSig = 0., valBkg = 0.;
1337  // If there's more than one shift, we should add them in quadrature
1338  if (fShifts.size() > 1) {
1339  for (auto& shift : fShifts) {
1340  double binShiftSig = shift.second->GetShiftedSpectrum(sigma).first->GetBinContent(bin) -
1341  fSpectra[0].first->GetBinContent(bin);
1342  double binShiftBkg = shift.second->GetShiftedSpectrum(sigma).second->GetBinContent(bin) -
1343  fSpectra[0].second->GetBinContent(bin);
1344  valSig += pow(binShiftSig,2);
1345  valBkg += pow(binShiftBkg,2);
1346  }
1347  double sign = sigma > 0 ? 1.0 : -1.0;
1348  shiftSig->SetBinContent(bin, fSpectra[0].first->GetBinContent(bin) + (sign * sqrt(valSig)));
1349  shiftBkg->SetBinContent(bin, fSpectra[0].second->GetBinContent(bin) + (sign * sqrt(valBkg)));
1350  }
1351  // If there's only one shift, just use the raw value
1352  else {
1353  valSig = fShifts.begin()->second->GetShiftedSpectrum(sigma).first->GetBinContent(bin)
1354  - fSpectra[0].first->GetBinContent(bin);
1355  valBkg = fShifts.begin()->second->GetShiftedSpectrum(sigma).second->GetBinContent(bin)
1356  - fSpectra[0].second->GetBinContent(bin);
1357  shiftSig->SetBinContent(bin, fShifts.begin()->second->GetShiftedSpectrum(sigma).first->GetBinContent(bin));
1358  shiftBkg->SetBinContent(bin, fShifts.begin()->second->GetShiftedSpectrum(sigma).second->GetBinContent(bin));
1359  }
1360  totShiftSig += valSig;
1361  totShiftBkg += valBkg;
1362  }
1363  fTotalShifts[sigma] = make_pair(totShiftSig, totShiftBkg);
1364  fSpectra[sigma] = make_pair(shiftSig, shiftBkg);
1365  }
1366 
1367  } // function SystMaker::ProcessSyst
1368 
1369  // -------------------------------------------------------------------------------------
1370  void SystMaker::SetAxes(const HistAxis* axis, const HistAxis* numuAxis) {
1371  fAxis = axis;
1372  fNuMuAxis = axis;
1373  }
1374 
1375  // -------------------------------------------------------------------------------------
1376  void SystMaker::SetCuts(const Cut* FDCut, const Cut* NDCut, const Cut* numuCut) {
1377  fNDCut = NDCut;
1378  fFDCut = FDCut;
1379  fNuMuCut = numuCut;
1380  }
1381 
1382  // -------------------------------------------------------------------------------------
1384  fLoaders = loaders;
1385  }
1386 
1387  // -------------------------------------------------------------------------------------
1389  fWeight = weight;
1390  }
1391 
1392  // -------------------------------------------------------------------------------------
1393  /// Function to make ISyst for this SystMaker
1395 
1396  if (fSpectra.size() == 0) {
1397  cout << "There are no spectra owned by this SystMaker!"
1398  << " Did you forget to call Go()?" << endl;
1399  return nullptr;
1400  }
1401 
1402  // Make ratios from shifted spectra
1403  map<int, pair<TH1*, TH1*>> shifts;
1404  for (const auto& sigma : fSpectra) {
1405  TH1* sigRatio = (TH1*)sigma.second.first->Clone(Form("SigRatioSigma%d", sigma.first));
1406  TH1* bkgRatio = (TH1*)sigma.second.second->Clone(Form("BgRatioSigma%d", sigma.first));
1407  //sigRatio->Divide(fSpectra.at(0).first);
1408  //bkgRatio->Divide(fSpectra.at(0).second);
1409  // BEGIN DEBUG CODE
1410  for (int i = 1; i <= sigRatio->GetNbinsX(); ++i) {
1411  double eps = 1e-3;
1412  double num = sigRatio->GetBinContent(i);
1413  double den = fSpectra.at(0).first->GetBinContent(i);
1414  double val = (num+eps)/(den+eps);
1415  // if (val > 2) {
1416  // cout << "Signal bin " << i << ", dividing thru "
1417  // << num << " by " << den << " to get " << val
1418  // << " - with epsilon it's " << valEps << endl;
1419  // }
1420  sigRatio->SetBinContent(i, val);
1421  num = bkgRatio->GetBinContent(i);
1422  den = fSpectra.at(0).second->GetBinContent(i);
1423  val = (num+eps)/(den+eps);
1424  bkgRatio->SetBinContent(i, val);
1425  // if (val > 2) {
1426  // cout << "Background bin " << i << ", dividing thru "
1427  // << num << " by " << den << " to get " << val
1428  // << " - with epsilon it's " << valEps << endl;
1429  // }
1430  }
1431  // END DEBUG CODE
1432  shifts[sigma.first] = make_pair(sigRatio, bkgRatio);
1433  }
1434 
1435  // Return new ISyst
1436  string shortName = fSampleName.empty() ? fName : fSampleName + "_" + fName;
1437  return new NuISyst(shortName, fLabel, fSampleName, shifts);
1438 
1439  } // function SystMaker::MakeISyst
1440 
1441  // -------------------------------------------------------------------------------------
1442  /// SaveTo implementation for SystMaker
1443  void SystMaker::SaveTo(TDirectory* dir) const {
1444 
1445  TDirectory* tmp = gDirectory;
1446 
1447  dir->cd();
1448  TObjString("SystMaker").Write("type");
1449 
1450  TObjString(fName.c_str()).Write("name");
1451  TObjString(fLabel.c_str()).Write("label");
1452  TObjString(fSampleName.c_str()).Write("samplename");
1453 
1454  fNominalPred->SaveTo(dir, "nompred");
1455 
1456  const char* shifts_dir = "Shifts";
1457  dir->mkdir("Shifts");
1458  for (const auto& shift : fShifts) {
1459  const char* shiftSubdir = Form("%s/%s", shifts_dir, shift.first.c_str());
1460  dir->mkdir(shiftSubdir);
1461  shift.second->SaveTo(dir->GetDirectory(shiftSubdir));
1462  }
1463 
1464  tmp->cd();
1465 
1466  } // function SystMaker::SaveTo
1467 
1468  // -------------------------------------------------------------------------------------
1469  /// LoadFrom implementation for SystMaker
1470  unique_ptr<SystMaker> SystMaker::LoadFrom(TDirectory* dir) {
1471 
1472  DontAddDirectory guard;
1473 
1474  TObjString* tag = (TObjString*)dir->Get("type");
1475  assert(tag);
1476  assert(tag->GetString() == "SystMaker");
1477 
1478  TObjString* name = (TObjString*)dir->Get("name");
1479  assert(name);
1480 
1481  TObjString* label = (TObjString*)dir->Get("label");
1482  assert(label);
1483 
1484  TObjString* sampleName = (TObjString*)dir->Get("samplename");
1485  assert(sampleName);
1486 
1487  unique_ptr<SystMaker> systMaker
1488  = make_unique<SystMaker>(name->GetString().Data(),
1489  label->GetString().Data(), sampleName->GetString().Data());
1490 
1491  // Load prediction
1492  tag = (TObjString*)dir->Get("nompred/type");
1493  assert(tag);
1494  if (tag->GetString() == "NDPredictionSterile")
1495  systMaker->fNominalPred = NDPredictionSterile::LoadFrom(dir, "nompred").release();
1496  else if (tag->GetString() == "FDPredictionSterile")
1497  systMaker->fNominalPred = FDPredictionSterile::LoadFrom(dir, "nompred").release();
1498  else if (tag->GetString() == "PredictionSterile")
1499  systMaker->fNominalPred = PredictionSterile::LoadFrom(dir, "nompred").release();
1500  else assert(false && "Prediction type not recognised!");
1501 
1502  // Load shift makers
1503  map<string, SystMakerShift*> shifts;
1504 
1505  TIter shift_next(dir->GetDirectory("Shifts")->GetListOfKeys());
1506  TKey* shift_key;
1507  while ((shift_key = (TKey*)shift_next())) {
1508  SystMakerShift* shift = SystMakerShift::LoadFrom(dir->GetDirectory(Form("Shifts/%s",
1509  shift_key->GetName()))).release();
1510  shifts[shift->GetName()] = shift;
1511  }
1512  systMaker->fShifts = shifts;
1513 
1514  return move(systMaker);
1515 
1516  } // function SystMaker::LoadFrom
1517 
1518  // -------------------------------------------------------------------------------------
1519  /// SystematicsMaker constructor
1521  fName(name),
1522  fPredType(predType),
1523  fAxis(nullptr),
1524  fNuMuAxis(nullptr),
1525  fNDCut(nullptr),
1526  fFDCut(nullptr),
1527  fNuMuCut(nullptr),
1528  fLoaders(nullptr),
1529  fLightLoaders(nullptr),
1530  fWeight(nullptr),
1531  fSample(NullSample)
1532  {
1533  vector<SystPredType> supportedPredTypes = {SystPredType::kExtrap, SystPredType::kFD, SystPredType::kND};
1534  if (find(supportedPredTypes.begin(), supportedPredTypes.end(), predType) == supportedPredTypes.end()) {
1535  cout << "SystematicsMaker error: Unsupported prediction type. Aborting..." << endl;
1536  abort();
1537  }
1538  } // SystematicsMaker constructor
1539 
1540  // -------------------------------------------------------------------------------------
1541  /// SystematicsMaker destructor
1543  for (auto systMaker : fSystMakers) {
1544  if (systMaker.second) delete systMaker.second;
1545  }
1546  } // function SystematicsMaker::~SystematicsMaker
1547 
1548  // -------------------------------------------------------------------------------------
1549  /// Get the names of all attached systematics
1551 
1552  vector<string> systNames;
1553  for (const auto& syst : fSystMakers)
1554  systNames.push_back(syst.first);
1555  return systNames;
1556 
1557  } // function SystematicsMaker::GetSystematicsNames
1558 
1559  // -------------------------------------------------------------------------------------
1560  /// Make ISyst objects for each attached systematic and return them
1561  vector<NuISyst*> SystematicsMaker::MakeISysts() {
1562 
1563  vector<NuISyst*> systs;
1564  for (const auto& maker : fSystMakers) {
1565  NuISyst* syst = maker.second->MakeISyst(fSample);
1566  systs.push_back(syst);
1567  }
1568  return systs;
1569 
1570  } // function SystematicsMaker::MakeISysts
1571 
1572  // -------------------------------------------------------------------------------------
1574 
1575  syst->SetAxes(fAxis, fNuMuAxis);
1576  syst->SetCuts(fFDCut, fNDCut, fNuMuCut);
1577  syst->SetNominalWeight(fWeight);
1578  if (syst->GetName().find("LightLevel") != string::npos)
1580  else
1581  syst->SetNominalLoaders(fLoaders);
1582  syst->MakePredictions(fPredType);
1583  fSystMakers[syst->GetName()] = syst;
1584  }
1585 
1586  // -------------------------------------------------------------------------------------
1587  void SystematicsMaker::PrintSysts(ostream& os) {
1588  for (const auto& syst : fSystMakers)
1589  syst.second->PrintSyst(os);
1590  return;
1591  }
1592 
1593  // -------------------------------------------------------------------------------------
1595  for (const auto& syst : fSystMakers)
1596  syst.second->DrawSyst(outDir);
1597  return;
1598  }
1599 
1600  // -------------------------------------------------------------------------------------
1601  void SystematicsMaker::Go(Sample sample, IOscCalc* calc) {
1602 
1603  fSample = sample;
1604  for (const auto& syst : fSystMakers) {
1605  syst.second->ProcessSyst(fSample, calc);
1606  fSystMakers[syst.first] = syst.second;
1607  }
1608  }
1609 
1610  // -------------------------------------------------------------------------------------
1611  void SystematicsMaker::SetAxes(const HistAxis* axis, const HistAxis* numuAxis) {
1612  fAxis = axis;
1613  fNuMuAxis = numuAxis;
1614  }
1615 
1616  // -------------------------------------------------------------------------------------
1617  void SystematicsMaker::SetCuts(const Cut* FDCut, const Cut* NDCut,
1618  const Cut* numuCut) {
1619  fNDCut = NDCut;
1620  fFDCut = FDCut;
1621  fNuMuCut = numuCut;
1622  }
1623 
1624  // -------------------------------------------------------------------------------------
1626  fLoaders = loaders;
1627  fLightLoaders = loaders_light;
1628  }
1629 
1630  // -------------------------------------------------------------------------------------
1632  fWeight = weight;
1633  }
1634 
1635  // -------------------------------------------------------------------------------------
1636  /// function to add another SystematicsMaker to this one
1638 
1639  map<string, SystMaker*> systMakers;
1640  for (const auto& systMaker : systMakers) {
1641  if (fSystMakers.count(systMaker.first) and !overwrite)
1642  continue;
1643  fSystMakers[systMaker.first] = systMaker.second;
1644  }
1645  } // function SystematicsMaker::Add
1646 
1647  // -------------------------------------------------------------------------------------
1648  void SystematicsMaker::SaveTo(TDirectory* dir, const string& name) const {
1649 
1650  TDirectory* tmp = gDirectory;
1651 
1652  dir = dir->mkdir(name.c_str()); // switch to subdir
1653  dir->cd();
1654 
1655  TObjString("SystematicsMaker").Write("type");
1656 
1657  TObjString(fName.c_str()).Write("name");
1658 
1659  int predType = -1;
1660  switch (fPredType) {
1661  case SystPredType::kExtrap:
1662  predType = 0;
1663  break;
1664  case SystPredType::kFD:
1665  predType = 1;
1666  break;
1667  case SystPredType::kND:
1668  predType = 2;
1669  break;
1670  default:
1671  break;
1672  }
1673  TVectorD predtype(1);
1674  predtype[0] = predType;
1675  predtype.Write("predtype");
1676 
1677  TDirectory* systMakersDir = dir->mkdir("SystMakers");
1678  for (const auto& systMaker : fSystMakers) {
1679  systMaker.second->SaveTo(systMakersDir->mkdir(systMaker.first.c_str()));
1680  }
1681 
1682  tmp->cd();
1683 
1684  } // function SystematicsMaker::SaveTo
1685 
1686  //----------------------------------------------------------------------
1687  unique_ptr<SystematicsMaker> SystematicsMaker::LoadFrom(
1688  TDirectory* dir, const string& name) {
1689 
1690  dir = dir->GetDirectory(name.c_str()); // switch to subdir
1691  assert(dir);
1692 
1693  DontAddDirectory guard;
1694 
1695  TObjString* tag = (TObjString*)dir->Get("type");
1696  assert(tag);
1697  assert(tag->GetString() == "SystematicsMaker");
1698 
1699  TObjString* tsName = (TObjString*)dir->Get("name");
1700  assert(tsName);
1701  string sName(tsName->GetString().Data());
1702 
1703  const TVectorD predtype = *(TVectorD*)dir->Get("predtype");
1705  switch (static_cast<int>(std::round(predtype[0]))) {
1706  case 0:
1707  predType = SystPredType::kExtrap;
1708  break;
1709  case 1:
1710  predType = SystPredType::kFD;
1711  break;
1712  case 2:
1713  predType = SystPredType::kND;
1714  break;
1715  }
1716 
1717  unique_ptr<SystematicsMaker> systsMaker = make_unique<SystematicsMaker>(sName, predType);
1718 
1719  map<string, SystMaker*> systMakers;
1720  TIter next(dir->GetDirectory("SystMakers")->GetListOfKeys());
1721  TKey* key;
1722  while ((key = (TKey*)next())) {
1723  SystMaker* systMaker
1724  = SystMaker::LoadFrom((TDirectory*)key->ReadObj()).release();
1725  systMakers[systMaker->GetName()] = systMaker;
1726  }
1727 
1728  systsMaker->fSystMakers = systMakers;
1729 
1730  return move(systsMaker);
1731 
1732  } // function SystematicsMaker::LoadFrom
1733 
1734 } // namespace ana
void SetAxes(const HistAxis *axis, const HistAxis *numuAxis=nullptr)
Definition: SystMaker.cxx:1370
const XML_Char * name
Definition: expat.h:151
static SystMakerShiftSigma * LoadFrom(TDirectory *dir)
Definition: SystMaker.cxx:193
std::string WriteTable(std::string name)
Write a table with up & down shifts for this systematic.
Definition: SystMaker.cxx:1229
static std::unique_ptr< FDPredictionSterile > LoadFrom(TDirectory *dir, const std::string &name)
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::map< int, std::pair< double, double > > fTotalShifts
Definition: SystMaker.h:242
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
void SetNominalLoaders(Loaders *loaders)
Definition: SystMaker.cxx:1383
set< int >::iterator it
Antineutrinos-only.
Definition: IPrediction.h:50
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
std::map< int, const Var * > fWeights
Definition: SystMaker.h:172
std::map< std::string, SystMakerShift * > fShifts
Definition: SystMaker.h:240
static std::unique_ptr< SystMakerShift > LoadFrom(TDirectory *dir)
LoadFrom implementation for SystMakerShift.
Definition: SystMaker.cxx:746
void SaveTo(TDirectory *dir) const
Definition: SystMaker.cxx:157
(&#39; appearance&#39;)
Definition: IPrediction.h:18
const HistAxis * fNuMuAxis
Definition: SystMaker.h:244
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:176
std::unique_ptr< IPrediction > Generate(Loaders &loaders, const SystShifts &shiftMC=kNoShift) const override
virtual void MakePredictions(SystPredType predType, const HistAxis *axis, const HistAxis *numuAxis, const Cut *fd_cut, const Cut *nd_cut, const Cut *numu_cut, const Var *weight, Loaders *loaders)
Base class implementation of MakePredictions, which just throws an error.
Definition: SystMaker.cxx:708
std::unique_ptr< IPrediction > Generate(Loaders &loaders, const SystShifts &shiftMC=kNoShift) const override
const Var weight
void PrintSysts(std::ostream &os=std::cout)
Definition: SystMaker.cxx:1587
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
IPrediction * fNominalPred
Definition: SystMaker.h:249
SystMakerWeightShift(std::string name, std::string label)
SystMakerWeightShift constructor.
Definition: SystMaker.cxx:848
Float_t den
Definition: plot.C:36
std::pair< TH1 *, TH1 * > fShiftedRatio
Definition: SystMaker.h:81
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
T sqrt(T number)
Definition: d0nt_math.hpp:156
std::pair< TH1 *, TH1 * > GetShiftedSpectrum(int sigma)
Definition: SystMaker.cxx:676
std::vector< std::string > GetSystematicsNames()
Get the names of all attached systematics.
Definition: SystMaker.cxx:1550
constexpr T pow(T x)
Definition: pow.h:72
std::vector< int > GetSigmas()
Get the sigmas for a given systematic shift.
Definition: SystMaker.cxx:982
SystMaker(std::string name, std::string label, std::string sample_name="")
SystMaker constructor.
Definition: SystMaker.cxx:923
const Cut * fFDCut
Definition: SystMaker.h:245
void PrintSyst(std::ostream &os=std::cout)
Print the total shifts for each shift owned by this syst maker.
Definition: SystMaker.cxx:1001
void PrintShift(std::ostream &os=std::cout)
Definition: SystMaker.cxx:256
static std::unique_ptr< PredictionSterile > LoadFrom(TDirectory *dir, const std::string &name)
~SystMakerWeightShift()
SystMakerWeightShift constructor.
Definition: SystMaker.cxx:854
TCanvas * canv
OStream cerr
Definition: OStream.cxx:7
virtual ~SystMakerShift()
Default destructor.
Definition: SystMaker.cxx:240
void abs(TH1 *hist)
void AddSigma(Loaders *loaders, int sigma)
Add sigma for loader shift with the same ND and FD loader.
Definition: SystMaker.cxx:781
const HistAxis * fAxis
Definition: SystMaker.h:244
SystMakerLoaderShift(std::string name, std::string label)
SystMakerLoaderShift constructor.
Definition: SystMaker.cxx:775
Float_t tmp
Definition: plot.C:36
virtual void SaveTo(TDirectory *dir, const std::string &name) const
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:336
SystMakerShift * GetShift(std::string shift_name)
Get the shift for a given name.
Definition: SystMaker.cxx:951
~SystMaker()
SystMaker destructor.
Definition: SystMaker.cxx:940
std::string fSampleName
Definition: SystMaker.h:238
Generates extrapolated NC predictions using ProportionalDecomp.
osc::OscCalcDumb calc
TH1 * GetSpectrum(IPrediction *pred, bool signal, const Sample sample, IOscCalc *calc)
Definition: SystMaker.cxx:76
void Add(SystematicsMaker *systMaker, bool overwrite)
function to add another SystematicsMaker to this one
Definition: SystMaker.cxx:1637
_IOscCalc< double > IOscCalc
std::string outDir
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:39
const Var * fWeight
Definition: SystMaker.h:295
void AddShift(SystMakerShift *shift)
Add SystMakerShift to this syst maker.
Definition: SystMaker.cxx:1015
const char * label
std::map< int, SystMakerShiftSigma * > fSigmas
Definition: SystMaker.h:128
void AddSigma(SystShifts *syst, int sigma)
Add sigma for syst shift.
Definition: SystMaker.cxx:901
void SetCuts(const Cut *fd_cut, const Cut *nd_cut, const Cut *numu_cut=nullptr)
Definition: SystMaker.cxx:1376
std::string GetName()
Definition: SystMaker.h:99
std::string fName
Definition: SystMaker.h:286
void SetAxes(const HistAxis *axis, const HistAxis *numuAxis=nullptr)
Definition: SystMaker.cxx:1611
const int nbins
Definition: cellShifts.C:15
std::map< int, std::pair< Loaders *, Loaders * > > fLoaders
Definition: SystMaker.h:151
Charged-current interactions.
Definition: IPrediction.h:39
TCanvas * DrawBkg()
Make plots of all the shifted spectra and ratios for this systmaker.
Definition: SystMaker.cxx:1109
std::vector< NuISyst * > MakeISysts()
Make ISyst objects for each attached systematic and return them.
Definition: SystMaker.cxx:1561
void MakePredictions(SystPredType predType)
Make predictions for all the SystMakerShift objects owned by this syst maker.
Definition: SystMaker.cxx:1296
static std::unique_ptr< NDPredictionSterile > LoadFrom(TDirectory *dir, const std::string &name)
void SetNominalLoaders(Loaders *loader, Loaders *loader_light=nullptr)
Definition: SystMaker.cxx:1625
Loaders * fLightLoaders
Definition: SystMaker.h:294
SystPredType fPredType
Definition: SystMaker.h:287
static std::unique_ptr< SystMaker > LoadFrom(TDirectory *dir)
LoadFrom implementation for SystMaker.
Definition: SystMaker.cxx:1470
void DrawSyst(TDirectory *outDir)
Make plots of all the shifted spectra and ratios for this systmaker.
Definition: SystMaker.cxx:1195
virtual void MakePredictions(SystPredType predType, const HistAxis *axis, const HistAxis *numuAxis, const Cut *fd_cut, const Cut *nd_cut, const Cut *numu_cut, const Var *weight, Loaders *loaders) override
Make shifted prediction for systematic shift derived from alternate MC loader.
Definition: SystMaker.cxx:826
void AddSigma(const Var *weight, int sigma)
Add sigma for weight shift.
Definition: SystMaker.cxx:861
void SaveTo(TDirectory *dir) const
SaveTo implementation for SystMaker.
Definition: SystMaker.cxx:1443
void SetCuts(const Cut *fd_cut, const Cut *nd_cut, const Cut *numu_cut=nullptr)
Definition: SystMaker.cxx:1617
covmx::Sample fSample
Definition: SystMaker.h:296
const Cut * fFDCut
Definition: SystMaker.h:292
virtual void MakePredictions(SystPredType predType, const HistAxis *axis, const HistAxis *numuAxis, const Cut *fd_cut, const Cut *nd_cut, const Cut *numu_cut, const Var *weight, Loaders *loaders) override
Make shifted prediction for systematic shift derived from SystShift.
Definition: SystMaker.cxx:910
std::string GetName()
Definition: SystMaker.h:203
void DrawSysts(TDirectory *outDir)
Definition: SystMaker.cxx:1594
virtual void MakePredictions(SystPredType predType, const HistAxis *axis, const HistAxis *numuAxis, const Cut *fd_cut, const Cut *nd_cut, const Cut *numu_cut, const Var *weight, Loaders *loaders) override
Make shifted prediction for systematic shift derived from CAFAna weight.
Definition: SystMaker.cxx:882
void AddSystematic(SystMaker *syst)
Definition: SystMaker.cxx:1573
static std::unique_ptr< SystematicsMaker > LoadFrom(TDirectory *dir, const std::string &name)
Definition: SystMaker.cxx:1687
IPrediction * fPrediction
Definition: SystMaker.h:79
void AddOnOff(Loaders *loaders)
Add on/off for loader shift.
Definition: SystMaker.cxx:802
std::string WriteTable(std::pair< TH1 *, TH1 * > nominal, std::string name)
Definition: SystMaker.cxx:609
void PrintSigma(std::ostream &os=std::cout)
Print the sigmas for a given shift.
Definition: SystMaker.cxx:149
std::vector< int > GetSigmas()
Definition: SystMaker.cxx:248
const HistAxis * fNuMuAxis
Definition: SystMaker.h:291
void ProcessSyst(const covmx::Sample sample, osc::IOscCalc *calc)
Definition: SystMaker.cxx:1313
std::string fName
Definition: SystMaker.h:125
float bin[41]
Definition: plottest35.C:14
~SystematicsMaker()
SystematicsMaker destructor.
Definition: SystMaker.cxx:1542
std::pair< double, double > fShift
Definition: SystMaker.h:82
std::string fLabel
Definition: SystMaker.h:126
double sigma(TH1F *hist, double percentile)
Neutrinos-only.
Definition: IPrediction.h:49
const SystShifts kNoShift
Definition: SystShifts.cxx:22
OStream cout
Definition: OStream.cxx:6
(&#39; survival&#39;)
Definition: IPrediction.h:19
const Cut * fNDCut
Definition: SystMaker.h:245
std::map< std::string, SystMaker * > fSystMakers
Definition: SystMaker.h:289
void ProcessShift(std::pair< TH1 *, TH1 * > &nominal_spectrum, const covmx::Sample sample, osc::IOscCalc *calc)
Definition: SystMaker.cxx:686
std::string fName
Definition: SystMaker.h:236
TString MakeLatexCommandName(TString mystring)
std::pair< TH1 *, TH1 * > fShiftedSpectrum
Definition: SystMaker.h:80
ifstream in
Definition: comparison.C:7
const Var * fWeight
Definition: SystMaker.h:247
TDirectory * dir
Definition: macro.C:5
~SystMakerShiftSigma()
SystMakerShiftSigma destructor.
Definition: SystMaker.cxx:139
std::map< std::string, std::string > systNames
Definition: PlotUnfolding.C:32
TCanvas * DrawSig()
Make plots of all the shifted spectra and ratios for this systmaker.
Definition: SystMaker.cxx:1023
int num
Definition: f2_nu.C:119
TCanvas * DrawSig(TH1 *nominal)
Definition: SystMaker.cxx:273
void DrawShift(TDirectory *outDir, std::pair< TH1 *, TH1 * > &nominal)
Definition: SystMaker.cxx:445
std::vector< Loaders * > loaders
Definition: syst_header.h:386
void SaveTo(TDirectory *dir) const
SaveTo implementation for SystMakerShift.
Definition: SystMaker.cxx:722
TCanvas * DrawBkg(TH1 *nominal)
Definition: SystMaker.cxx:358
void Go(const covmx::Sample sample, osc::IOscCalc *calc)
Definition: SystMaker.cxx:1601
Neutral-current interactions.
Definition: IPrediction.h:40
void SetNominalWeight(const Var *weight)
Definition: SystMaker.cxx:1388
assert(nhit_max >=nhit_nbins)
string shortName
THUMBNAIL BLOCK: We need to make a thumbnail for each.
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
std::unique_ptr< IPrediction > Generate(Loaders &loaders, const SystShifts &shiftMC=kNoShift) const override
const Cut * fNuMuCut
Definition: SystMaker.h:245
IPrediction * GetPrediction(SystPredType predType, const HistAxis *axis, const HistAxis *numuAxis, const Cut *FDCut, const Cut *NDCut, const Cut *numuCut, const SystShifts *shiftData, const SystShifts *shiftMC, const Var *weight, Loaders *loaders, Loaders *loadersND)
Definition: SystMaker.cxx:42
Loaders * fLoaders
Definition: SystMaker.h:246
void SetNominalWeight(const Var *weight)
Definition: SystMaker.cxx:1631
std::map< int, std::pair< double, double > > GetTotalShift()
SystematicsMaker(std::string name, SystPredType predType)
SystematicsMaker constructor.
Definition: SystMaker.cxx:1520
std::string fLabel
Definition: SystMaker.h:237
void AddOnOff(const Var *weight)
Add on/off for weight shift.
Definition: SystMaker.cxx:870
SystMakerSystShift(std::string name, std::string label)
SystMakerSystShift constructor.
Definition: SystMaker.cxx:895
std::map< int, std::pair< TH1 *, TH1 * > > fSpectra
Definition: SystMaker.h:241
All neutrinos, any flavor.
Definition: IPrediction.h:26
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
Float_t e
Definition: plot.C:35
Generates Near Detector predictions.
NuISyst * MakeISyst(covmx::Sample sample)
Function to make ISyst for this SystMaker.
Definition: SystMaker.cxx:1394
void next()
Definition: show_event.C:84
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: SystMaker.cxx:1648
SystMakerShift(std::string name, std::string label)
Definition: SystMaker.cxx:233
const HistAxis * fAxis
Definition: SystMaker.h:291
def sign(x)
Definition: canMan.py:197
const Cut * fNDCut
Definition: SystMaker.h:292
const Cut * fNuMuCut
Definition: SystMaker.h:292
gm Write()
SystPredType
Definition: SystMaker.h:38
std::map< int, SystShifts * > fSysts
Definition: SystMaker.h:190
const covmx::Sample NullSample(covmx::Selection::kNoSel, covmx::Polarity::kNoPol, covmx::Detector::kNoDet)
enum BeamMode string