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