FDDataMCSystBandAna.C
Go to the documentation of this file.
1 // Macro to plot FD Data/MC spectra with systematic error band
2 
5 #include "CAFAna/Core/Loaders.h"
6 #include "CAFAna/Core/Spectrum.h"
7 #include "CAFAna/Cuts/Cuts.h"
8 #include "NuXAna/Cuts/NusCuts.h"
13 #include "NuXAna/Systs/NusSysts.h"
14 #include "CAFAna/Systs/Systs.h"
17 
18 #include "TCanvas.h"
19 #include "TFile.h"
20 #include "TGraphAsymmErrors.h"
21 #include "TH1.h"
22 #include "TLegend.h"
23 #include "TStyle.h"
24 
25 #include <iostream>
26 #include <string>
27 #include <vector>
28 
29 using namespace ana;
30 
32  const std::vector<const ISyst*>& systs,
34  TDirectory* out, FILE* text, strings strs,
35  Spectrum& sCos, const Spectrum& sOOT,
36  Spectrum* sData);
37 
38 void PlotSysts(Spectrum spectra[],
39  const std::vector<Spectrum>& upShifts,
40  const std::vector<Spectrum>& downShifts,
41  TDirectory* out, FILE* text, strings strs,
42  double mcscale, const Spectrum& sOOT,
43  Spectrum* sData);
44 
46  const ISyst* syst,
48  TDirectory* out, strings strs);
49 
50 void PlotSysts(Spectrum& nom, Spectrum& sp1, Spectrum& sm1,
51  TDirectory* out, strings strs, bool split = true);
52 
53 TH1* GetHist(strings strs, int sigma);
54 
56 {
57  TH1::AddDirectory(0);
58 
59  // Create vector of systematics
60  std::vector<const ISyst*> systs = getAllNusSysts();
61 
62  std::string folder = "/nova/ana/steriles/Ana01/";
63  std::string filenm = "FDDataMCSystBand";
64 
65  std::string loadLocation = folder + filenm + ".root";
66  std::string saveLocation = folder + filenm + "Ana.root";
67  std::string textLocation = folder + filenm + "Ana.txt";
68 
69  TFile* rootL = new TFile(loadLocation.c_str(), "READ");
70  TDirectory* tmpL = gDirectory;
71  TDirectory* loadDir = gDirectory;
72 
73  loadDir->cd((loadLocation + ":/pred").c_str());
75 
76  loadDir->cd((loadLocation + ":/sCosmic").c_str());
77  Spectrum sCosmic = *Spectrum::LoadFrom(gDirectory);
78 
79  loadDir->cd((loadLocation + ":/sCosOOT").c_str());
80  Spectrum sCosOOT = *Spectrum::LoadFrom(gDirectory);
81 
82  loadDir->cd((loadLocation + ":/sData").c_str());
83  Spectrum sData = *Spectrum::LoadFrom(gDirectory);
84 
85  tmpL->cd();
86  rootL->Close(); // Don't forget to close the file!
87 
88  TFile* rootF = new TFile(saveLocation.c_str(), "RECREATE");
89  FILE* textF;
90  textF = fopen(textLocation.c_str(), "w");
91  fprintf(textF, "FD Systematics\n\n");
92 
93  strings strs;
94  strs.fComponent = "CalE";
95  strs.fDet = "EX";
96  strs.fPOT = "";
97  strs.fSample = "";
98  strs.fSystType = "";
99  strs.fXLabel = "Calorimetric Energy (GeV)";
100 
101  strs.fSystS = "All";
102  strs.fSystL = "All";
103 
105  PlotSysts(&pred, systs, calc, rootF, textF, strs, sCosmic, sCosOOT, &sData);
106 
107  for(unsigned int i = 0; i < systs.size(); ++i)
108  {
109  strs.fSystS = systs[i]->ShortName();
110  strs.fSystL = systs[i]->LatexName();
111 
112  PlotSysts(&pred, {systs[i]}, calc, rootF, textF, strs, sCosmic, sCosOOT, &sData);
113  PlotSysts(&pred, systs[i], calc, rootF, strs);
114  }
115 
116  fclose(textF);
117  rootF->Close();
118 }
119 
120 //----------------------------------------------------------------------
122  const std::vector<const ISyst*>& systs,
124  TDirectory* out, FILE* text, strings strs,
125  Spectrum& sCos, const Spectrum& sOOT,
126  Spectrum* sData)
127 {
128  if(!sCos.POT()) {
130  }
131 
132  double mcscale = 0.;
133  mcscale += pred->GetPrediction(0)->Predict(calc).POT();
134  mcscale += pred->GetPrediction(1)->Predict(calc).POT();
135 
136  Spectrum spectraFD[MAXSPEC] = {
137  sCos + pred->Predict(calc),
141  sCos
142  };
143 
144  std::vector<Spectrum> ups, dns;
145 
146  for(const ISyst* syst: systs){
147  SystShifts shifts;
148  shifts.SetShift(syst, +1);
149  ups.push_back(pred->PredictSyst(calc, shifts));
150  shifts.SetShift(syst, -1);
151  dns.push_back(pred->PredictSyst(calc, shifts));
152  }
153 
154  PlotSysts(spectraFD, ups, dns, out, text, strs, mcscale, sOOT, sData);
155 }
156 
157 //----------------------------------------------------------------------
158 void PlotSysts(Spectrum spectra[],
159  const std::vector<Spectrum>& upShifts,
160  const std::vector<Spectrum>& downShifts,
161  TDirectory* out, FILE* text, strings strs,
162  double mcscale, const Spectrum& sOOT,
163  Spectrum* dat)
164 {
165  double scale = kSecondAnaPOT;
166 
167  // Get nominal histograms
168 
169  // Get each histogram from the Spectrum array
170  TH1* hmc = spectra[0].ToTH1(scale);
171  TH1* hNC = spectra[1].ToTH1(scale);
172  TH1* hNumu = spectra[2].ToTH1(scale);
173  TH1* hNue = spectra[3].ToTH1(scale);
174  TH1* hCos = spectra[4].ToTH1(scale);
175  TH1* hdt = dat->ToTH1(dat->POT());
176  TH1* hOOT = sOOT.ToTH1(sOOT.Livetime(), kLivetime);
177 
178  // Set up vectors of shifted histograms
179  std::vector<TH1*> ups, dns;
180  for(const auto& upShift:upShifts) ups.push_back(upShift.ToTH1(scale));
181  for(const auto& downShift:downShifts) dns.push_back(downShift.ToTH1(scale));
182 
183  // Construct histogram names and titles
184  std::string xlabel = strs.fXLabel;
185  std::string ylabel = "Events / 0.25 GeV";
186  std::string title = strs.fDet + " " + strs.fComponent +
187  " Error for " + strs.fSystL + " Systematic";
188 // std::string fullTitle = title + ";" + xlabel + ";" + ylabel;
189  std::string fullTitle = ";" + xlabel + ";" + ylabel;
190  std::string name = "c" + strs.fComponent + strs.fDet + strs.fSystS;
191  std::string ylabelRat = "Ratio to Nominal MC";
192  std::string fullTitleRat = ";" + xlabel + ";" + ylabelRat;
193 
194  fprintf(text, "%s\n", strs.fSystL.c_str());
195 
196  double tot_events = 0.;
197  double up_events = 0.;
198  double dn_events = 0.;
199 
200  // Variables for setting axis ranges
201  double maxval = 0.;
202  double maxvalRat = 0.;
203  double minvalRat = 1.;
204 
205  // Set up histograms for ratio canvas
206  TH1* hDtRat = (TH1*)hdt->Clone();
207  hDtRat->Divide(hmc);
208  TH1* hMCRat = (TH1*)hmc->Clone();
209 
210  // Set up objects for containing the MC errors
211  TGraphAsymmErrors* gmc = new TGraphAsymmErrors;
212  TGraphAsymmErrors* gMCRat = new TGraphAsymmErrors;
213  TGraphAsymmErrors* gFullErr = new TGraphAsymmErrors;
214 
215  // Loop over bins, don't include under/overflow bins
216  for(int i_bin = 1, n_bin = hmc->GetNbinsX(); i_bin <= n_bin; ++i_bin) {
217  double valmc = hmc->GetBinContent(i_bin);
218  double valcs = hCos->GetBinContent(i_bin);
219  double valdt = hdt->GetBinContent(i_bin);
220  tot_events += valmc;
221 
222  const double w = hmc->GetXaxis()->GetBinWidth(i_bin);
223  double errUp = 0.;
224  double errDn = 0.;
225  double errSt = 0.;
226 
227  // Set the bin values
228  hMCRat->SetBinContent(i_bin, 1.);
229  gmc ->SetPoint(i_bin, hmc->GetXaxis()->GetBinCenter(i_bin), valmc);
230  gFullErr->SetPoint(i_bin, hmc->GetXaxis()->GetBinCenter(i_bin), valmc);
231  gMCRat ->SetPoint(i_bin, hmc->GetXaxis()->GetBinCenter(i_bin), 1.);
232 
233  for(unsigned int i_syst = 0; i_syst < ups.size(); ++i_syst) {
234  double hi = ups[i_syst]->GetBinContent(i_bin) + valcs - valmc;
235  double lo = dns[i_syst]->GetBinContent(i_bin) + valcs - valmc; // This should be negative!
236 
237  // Correct values under a few "error" conditions
238  if(hi < 0. && lo > 0.) {
239  std::swap(lo, hi);
240  }
241  if(hi < 0.) {
242  lo = std::min(lo, hi);
243  hi = 0.;
244  }
245  if(lo > 0.) {
246  hi = std::max(lo, hi);
247  lo = 0.;
248  }
249 
250  errUp += hi*hi;
251  errDn += lo*lo;
252  } // End loop over systs
253 
254  errUp = sqrt(errUp);
255  errDn = sqrt(errDn);
256 
257  errSt += sqrt(hNC ->GetBinContent(i_bin));
258  errSt += sqrt(hNue ->GetBinContent(i_bin));
259  errSt += sqrt(hNumu->GetBinContent(i_bin));
260  errSt *= sqrt(scale/mcscale);
261  if(hCos->GetBinContent(i_bin) > 0) {
262  errSt += hCos->GetBinContent(i_bin)/sqrt(hOOT->GetBinContent(i_bin));
263  }
264 
265  up_events += errUp;
266  dn_events += errDn;
267 
268  // Calculate a few quantities and protect against Inf/NaN
269  double errUpPct = ( (valmc > 0.) ? errUp/valmc : 0. );
270  double errDnPct = ( (valmc > 0.) ? errDn/valmc : 0. );
271  double errData = sqrt(valdt);
272  double errDataPct = ( (valdt > 0.) ? errData/valdt : 0. );
273  double dtrat = ( (valmc > 0.) ? valdt/valmc : 1. );
274 
275  // Set bin errors
276  hdt ->SetBinError(i_bin, errData);
277  hDtRat->SetBinError(i_bin, errDataPct);
278 
279  gmc ->SetPointError(i_bin, w/2, w/2, errDn, errUp);
280  gMCRat ->SetPointError(i_bin, w/2, w/2, errDnPct, errUpPct);
281  gFullErr->SetPointError(i_bin, w/2, w/2, errDn + errSt, errUp + errSt);
282 
283  // Calculate maximum values for axis ranges
284  maxval = std::max(maxval, valmc + errUp);
285  maxval = std::max(maxval, valdt + errData);
286 
287  maxvalRat = std::max(maxvalRat, 1. + errUpPct);
288  maxvalRat = std::max(maxvalRat, dtrat + errDataPct);
289 
290  minvalRat = std::min(minvalRat, 1. - errDn/valmc);
291  minvalRat = std::min(minvalRat, dtrat - errDataPct);
292  } // End loop over bins
293 
294  minvalRat = 1. - std::abs(1. - minvalRat);
295  maxvalRat = 1. + std::abs(1. - maxvalRat);
296 
297  SetHistOptions(hmc, maxval, fullTitle, -404, kTotalMCColor, false);
298  gmc->SetFillColor(kTotalMCErrorBandColor);
299  gmc->SetLineColor(kTotalMCColor);
300  gmc->SetLineWidth(2);
301  gFullErr->SetFillColor(kTotalMCErrorBandColor);
302  gFullErr->SetLineColor(kTotalMCColor);
303  gFullErr->SetLineWidth(2);
304 
305  CenterTitles(hMCRat);
306  hMCRat->SetLineColor(kTotalMCColor);
307  hMCRat->SetMaximum(1.05*maxvalRat);
308  hMCRat->SetMinimum(0.95*minvalRat);
309  hMCRat->SetTitle(fullTitleRat.c_str());
310  gMCRat->SetLineColor(kTotalMCColor);
311  gMCRat->SetFillColor(kTotalMCErrorBandColor);
312  gMCRat->SetLineWidth(2);
313 
314  SetHistOptions(hdt, maxval, fullTitle, -404, kBlack, false);
315  hdt->SetMarkerStyle(kFullCircle);
316 
317  const Color_t kCosmicBackgroundColor = kOrange+7;
318  SetHistOptions(hNC, maxval, fullTitle, -404, kAzure+2, false);
319  SetHistOptions(hNue, maxval, fullTitle, -404, kBeamNueBackgroundColor, false);
320  SetHistOptions(hNumu, maxval, fullTitle, -404, kNumuBackgroundColor, false);
321  SetHistOptions(hCos, maxval, fullTitle, -404, kCosmicBackgroundColor, false);
322 
323  CenterTitles(hDtRat);
324  hDtRat->SetLineColor(kBlack);
325  hDtRat->SetMaximum(1.05*maxvalRat);
326  hDtRat->SetMinimum(0.95*minvalRat);
327  hDtRat->SetTitle(fullTitleRat.c_str());
328 
329  double xL = 0.6, xR = 0.85;
330  double yB = 0.6, yT = 0.85;
331  TLegend* leg = new TLegend(xL, yB, xR, yT);
332  SetLegendOptions(leg);
333  leg->AddEntry(hdt, "FD Data", "lep");
334  leg->AddEntry(gmc, "Total MC Prediction", "lf");
335  leg->SetY1(leg->GetY2() - leg->GetNRows()*0.05);
336 
337  // Plot Data/MC, just spectrum
338  TCanvas* cNoRat = new TCanvas(name.c_str(), title.c_str(), 800, 500);
339  gPad->SetFillStyle(0);
340  hmc->Draw("hist");
341  gmc->Draw("e2 same");
342  hmc->Draw("hist same");
343  hdt->Draw("same");
344  leg->Draw();
345  DrawLatexLines(xL + 0.01, leg->GetY1() - 0.045, "6.05");
346  cNoRat->RedrawAxis();
347 
348  Preliminary();
349  out->WriteTObject(cNoRat);
350 
351  // Plot Data/MC, just spectrum, add individual components
352  TLegend* legAll = new TLegend(xL, yB, xR, yT);
353  SetLegendOptions(legAll);
354  legAll->AddEntry(hdt, "FD Data", "lep");
355  legAll->AddEntry(gmc, "Total Prediction", "lf");
356  legAll->AddEntry(hNC, "NC Prediction", "l");
357  legAll->AddEntry(hNue, "#nu_{e} CC Background", "l");
358  legAll->AddEntry(hNumu, "#nu_{#mu} CC Background", "l");
359  legAll->AddEntry(hCos, "Cosmic Background", "l");
360  legAll->SetY1(legAll->GetY2() - legAll->GetNRows()*0.05);
361 
362  TCanvas* cAllMC = new TCanvas((name + "Comp").c_str(), (title + "Comp").c_str(), 800, 500);
363  gPad->SetFillStyle(0);
364  hmc ->Draw("hist");
365  gmc ->Draw("e2 same");
366  hmc ->Draw("hist same");
367  hNC ->Draw("hist same");
368  hCos ->Draw("hist same");
369  hNue ->Draw("hist same");
370  hNumu->Draw("hist same");
371  hdt ->Draw("same");
372  legAll->Draw();
373  DrawLatexLines(xL + 0.01, legAll->GetY1() - 0.045, "6.05");
374  cAllMC->RedrawAxis();
375 
376  Preliminary();
377  out->WriteTObject(cAllMC);
378 
379  // Plot Data/MC, just spectrum, add individual components, Stat and Systs
380  // There's probably a better way to do this, but this hack works.
381  TLegend* leg1 = new TLegend(xL, yB, xR, yT);
382  SetLegendOptions(leg1);
383  leg1->AddEntry(hdt, "FD Data", "lep");
384  leg1->SetY1(leg1->GetY2() - leg1->GetNRows()*0.05);
385 
386  TLegend* leg2 = new TLegend(xL, yB, xR, leg1->GetY1() - 0.015);
387  SetLegendOptions(leg2);
388  leg2->AddEntry(gFullErr, "#splitline{Total Prediction}{Stat. and Syst. Uncert.}", "lf");
389  leg2->SetY1(leg2->GetY2() - leg2->GetNRows()*0.05);
390 
391  TLegend* leg3 = new TLegend(xL, yB, xR, leg2->GetY1() - 0.0175);
392  SetLegendOptions(leg3);
393  leg3->AddEntry(hNC, "NC Prediction", "l");
394  leg3->AddEntry(hNue, "#nu_{e} CC Background", "l");
395  leg3->AddEntry(hNumu, "#nu_{#mu} CC Background", "l");
396  leg3->AddEntry(hCos, "Cosmic Background", "l");
397  leg3->SetY1(leg3->GetY2() - leg3->GetNRows()*0.05);
398 
399  TCanvas* cError = new TCanvas((name + "StSy").c_str(), (title + "StSy").c_str(), 800, 500);
400  gPad->SetFillStyle(0);
401  hmc->Draw("hist");
402  gFullErr->Draw("e2 same");
403  hmc->Draw("hist same");
404  hNC ->Draw("hist same");
405  hCos ->Draw("hist same");
406  hNue ->Draw("hist same");
407  hNumu->Draw("hist same");
408  hdt->Draw("same");
409  leg1->Draw();
410  leg2->Draw();
411  leg3->Draw();
412  DrawLatexLines(xL + 0.01, leg3->GetY1() - 0.045, "6.05");
413  cError->RedrawAxis();
414 
415  Preliminary();
416  out->WriteTObject(cError);
417 
418  TCanvas* c = new TCanvas((name + "Rat").c_str(), (title + " Ratio").c_str(), 800, 800);
419  gPad->SetFillStyle(0);
420  TPad* pSpecs = new TPad("pSpecs", "", 0., 0.375, 1., 1.);
421  TPad* pRatio = new TPad("pRatio", "", 0., 0., 1., 0.375);
422  pSpecs->Draw();
423  pRatio->Draw();
424 
425  pSpecs->cd();
426  hmc->Draw("hist");
427  gmc->Draw("e2 same");
428  hmc->Draw("hist same");
429  hdt->Draw("same");
430  leg->Draw();
431  pSpecs->RedrawAxis();
432  Preliminary();
433 
434  pRatio->cd();
435  hMCRat->Draw("hist");
436  gMCRat->Draw("e2 same");
437  hMCRat->Draw("hist same");
438  hDtRat->Draw("same");
439  pRatio->RedrawAxis();
440 
441  c->Update();
442  out->WriteTObject(c);
443 
444  fprintf(text, "Number of events: %6.2f, Error up: %6.2f, Error down: %6.2f\n",
445  tot_events, up_events, dn_events);
446  fprintf(text, "Percentage up: %4.2f%%, Percentage down: %4.2f%%",
447  100.*up_events/tot_events, 100.*dn_events/tot_events);
448  fprintf(text, "\n\n");
449 
450  return;
451 }
452 
453 //----------------------------------------------------------------------
455  const ISyst* syst,
457  TDirectory* out, strings strs)
458 {
461 
462  std::vector<Spectrum> ups, dns;
463 
464  SystShifts shifts;
465  shifts.SetShift(syst, +1);
468  shifts.SetShift(syst, -1);
471 
472  strs.fComponent = "NC";
473  PlotSysts(nomNC, upNC, dnNC, out, strs);
474  strs.fComponent = "BG";
475  PlotSysts(nomCC, upCC, dnCC, out, strs);
476 }
477 
478 //----------------------------------------------------------------------
479 void PlotSysts(Spectrum& nom, Spectrum& sp1, Spectrum& sm1,
480  TDirectory* out, strings strs, bool split)
481 {
482  // Scale everything to data POT
483  double scale = kSecondAnaPOT;
484 
485  TH1* h = nom.ToTH1(scale);
486  TH1* hp1 = sp1.ToTH1(scale);
487  TH1* hm1 = sm1.ToTH1(scale);
488 
489  // Construct histogram names and titles
490  std::string xlabel = strs.fXLabel;
491  std::string ylabel = "Events / 6.05 #times 10^{20} POT";
492  std::string title = strs.fDet + " " + strs.fComponent +
493  " Error for " + strs.fSystL + " Systematic";
494  std::string fullTitle = title + ";" + xlabel + ";" + ylabel;
495  std::string name = "c" + strs.fComponent + strs.fDet + strs.fSystS;
496  std::string ylabelRat = "Shifted / Nominal";
497  std::string fullTitleRat = ";" + xlabel + ";" + ylabelRat;
498 
499  TH1* hRat = (TH1*)h->Clone();
500  TH1* hp1Rat = (TH1*)hp1->Clone();
501  TH1* hm1Rat = (TH1*)hm1->Clone();
502 
503  int nBins = h->GetNbinsX();
504  for(int i = 0; i <= nBins+1; ++i) {
505  hRat->SetBinContent(i, 1.);
506  }
507  hp1Rat->Divide(h);
508  hm1Rat->Divide(h);
509 
510  double maxval = h->GetBinContent(h->GetMaximumBin());
511  maxval = std::max(maxval, hp1->GetBinContent(hp1->GetMaximumBin()));
512  maxval = std::max(maxval, hm1->GetBinContent(hm1->GetMaximumBin()));
513 
514  double minval = h->GetBinContent(h->GetMinimumBin());
515  minval = std::min(minval, hp1->GetBinContent(hp1->GetMinimumBin()));
516  minval = std::min(minval, hm1->GetBinContent(hm1->GetMinimumBin()));
517 
518  if(maxval < minval) { std::swap(maxval, minval); }
519 
520  double maxvalRat = hRat->GetBinContent(hRat->GetMaximumBin());
521  maxvalRat = std::max(maxvalRat, hp1Rat->GetBinContent(hp1Rat->GetMaximumBin()));
522  maxvalRat = std::max(maxvalRat, hm1Rat->GetBinContent(hm1Rat->GetMaximumBin()));
523 
524  double minvalRat = hRat->GetBinContent(hRat->GetMinimumBin());
525  for(int i = 1; i <= nBins; ++i) {
526  if(hp1Rat->GetBinContent(i) > 0.) {
527  minvalRat = std::min(minvalRat, hp1Rat->GetBinContent(i));
528  }
529  if(hm1Rat->GetBinContent(i) > 0.) {
530  minvalRat = std::min(minvalRat, hm1Rat->GetBinContent(i));
531  }
532  }
533 
534  if(std::abs(maxvalRat - 1.) > std::abs(1. - minvalRat)) {
535  minvalRat = 1 - std::abs(maxvalRat - 1.);
536  maxvalRat = 1 + std::abs(maxvalRat - 1.);
537  }
538  else {
539  minvalRat = 1 - std::abs(1. - minvalRat);
540  maxvalRat = 1 + std::abs(1. - minvalRat);
541  }
542 
543  CenterTitles(h);
544  h->SetLineColor(kBlack);
545  h->SetMaximum(1.1*maxval);
546  h->SetMinimum(0);
547  h->SetTitle(fullTitle.c_str());
548 
549  CenterTitles(hp1);
550  hp1->SetLineColor(kRed);
551  hp1->SetLineStyle(2);
552  hp1->SetMaximum(1.1*maxval);
553  hp1->SetMinimum(0);
554  hp1->SetTitle(fullTitle.c_str());
555 
556  CenterTitles(hm1);
557  hm1->SetLineColor(kBlue);
558  hm1->SetLineStyle(2);
559  hm1->SetMaximum(1.1*maxval);
560  hm1->SetMinimum(0);
561  hm1->SetTitle(fullTitle.c_str());
562 
563  CenterTitles(hRat);
564  hRat->SetLineColor(kBlack);
565  hRat->SetMaximum(1.1*maxvalRat);
566  hRat->SetMinimum(0.9*minvalRat);
567  hRat->SetTitle(fullTitleRat.c_str());
568 
569  CenterTitles(hp1Rat);
570  hp1Rat->SetLineColor(kRed);
571  hp1Rat->SetLineStyle(2);
572  hp1Rat->SetMaximum(1.1*maxvalRat);
573  hp1Rat->SetMinimum(0.9*minvalRat);
574  hp1Rat->SetTitle(fullTitleRat.c_str());
575 
576  CenterTitles(hm1Rat);
577  hm1Rat->SetLineColor(kBlue);
578  hm1Rat->SetLineStyle(2);
579  hm1Rat->SetMaximum(1.1*maxvalRat);
580  hm1Rat->SetMinimum(0.9*minvalRat);
581  hm1Rat->SetTitle(fullTitleRat.c_str());
582 
583  double xL = 0.6, xR = 0.85;
584  double yB = 0.6, yT = 0.85;
585  TLegend* leg = new TLegend(xL, yB, xR, yT);
586  leg->AddEntry(h, "Nominal", "l");
587  leg->AddEntry(hp1, "+1 Sigma", "l");
588  leg->AddEntry(hm1, "-1 Sigma", "l");
589 
590  TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 800);
591  gPad->SetFillStyle(0);
592  TPad* pSpecs = new TPad("pSpecs", "", 0., 0.375, 1., 1.);
593  TPad* pRatio = new TPad("pRatio", "", 0., 0., 1., 0.375);
594  pSpecs->Draw();
595  pRatio->Draw();
596 
597  pSpecs->cd();
598  h->Draw("hist");
599  hp1->Draw("hist same");
600  hm1->Draw("hist same");
601  leg->Draw();
602  Simulation();
603 
604  pRatio->cd();
605  hRat->Draw("hist");
606  hp1Rat->Draw("hist same");
607  hm1Rat->Draw("hist same");
608 
609  c->Update();
610  out->WriteTObject(c);
611 
612  if(split) {
613  TH1* hp1DRat = GetHist(strs, +1);
614  TH1* hm1DRat = GetHist(strs, -1);
615 
616  hp1DRat->Divide(hp1Rat);
617  hm1DRat->Divide(hm1Rat);
618 
619  CenterTitles(hp1DRat);
620  hp1DRat->SetLineColor(kBlack);
621  hp1DRat->SetMaximum(1.1*maxvalRat);
622  hp1DRat->SetMinimum(0.9*minvalRat);
623  hp1DRat->SetTitle(fullTitle.c_str());
624 
625  CenterTitles(hp1);
626  hm1DRat->SetLineColor(kRed);
627  hm1DRat->SetLineStyle(2);
628  hm1DRat->SetMaximum(1.1*maxvalRat);
629  hm1DRat->SetMinimum(0.9*minvalRat);
630  hm1DRat->SetTitle(fullTitle.c_str());
631 
632  TLegend* leg2 = new TLegend(xL, yB, xR, yT);
633  leg2->AddEntry(hp1DRat, "+1 Sigma Double Ratio", "l");
634  leg2->AddEntry(hm1DRat, "-1 Sigma Double Ratio", "l");
635 
636  std::string titleDR = strs.fDet + " " + strs.fComponent +
637  strs.fSystL + " Double Ratio";
638  std::string nameDR = name + "DR";
639  std::string ylabelDRat = "File Ratio / Macro Ratio";
640  std::string fullTitleDRat = ";" + xlabel + ";" + ylabelDRat + titleDR;
641 
642  TCanvas* cDR = new TCanvas(nameDR.c_str(), fullTitleDRat.c_str(), 800, 500);
643  gPad->SetFillStyle(0);
644  hp1DRat->Draw("hist");
645  hm1DRat->Draw("hist same");
646  leg2->Draw();
647  Simulation();
648 
649  out->WriteTObject(cDR);
650  }
651 
652  return;
653 }
654 
655 //----------------------------------------------------------------------
656 TH1* GetHist(strings strs, int sigma)
657 {
658  std::string filename = "$NOVA_ANA/steriles/Ana01/Systematics/NusSystsAna01.root";
659  TFile fin(filename.c_str(), "read");
660  if(fin.IsZombie()) {
661  std::cout << "Warning: couldn't open file." << std::endl;
662  abort();
663  }
664 
665  std::string syst = strs.fSystS;
666 
667  if(syst.compare("mcstat") && syst.compare("ndcont") && syst.compare("norm")) {
668  std::string strsigma = ((sigma == 1) ? "+1" : "-1");
669  std::string name = TString::Format("h%s_%s_%s_%s",
670  strs.fComponent.c_str(),
671  strs.fDet.c_str(),
672  strs.fSystS.c_str(),
673  strsigma.c_str()).Data();
674 
675  TH1* ret = (TH1*)fin.Get(name.c_str());
676 
677  if(!ret) {
678  std::cout << "Error: can't find necessary " << name << std::endl;
679  abort();
680  }
681 
682  return ret;
683  }
684 
685  std::string name = TString::Format("hNC_ND_Beam_0").Data();
686  TH1* ret = (TH1*)fin.Get(name.c_str());
687 
688  if(!ret) {
689  std::cout << "Error: can't find necessary " << name << std::endl;
690  abort();
691  }
692 
693  double shift = 1.;
694  double dsigma = (double)sigma;
695  std::string comp = strs.fComponent;
696  if(!comp.compare("NC") && !syst.compare("mcstat")) shift += dsigma*.020;
697  if(!comp.compare("BG") && !syst.compare("mcstat")) shift += dsigma*.048;
698  if(!comp.compare("NC") && !syst.compare("ndcont")) shift += dsigma*.010;
699  if(!comp.compare("BG") && !syst.compare("ndcont")) shift += dsigma*.006;
700  if(!comp.compare("NC") && !syst.compare("norm") ) shift += dsigma*.032;
701  if(!comp.compare("BG") && !syst.compare("norm") ) shift += dsigma*.031;
702 
703  for(int i = 0, n = ret->GetNbinsX(); i <= n+1; ++i) {
704  ret->SetBinContent(i, shift);
705  }
706 
707  return ret;
708 }
void Simulation()
Definition: tools.h:16
void SetHistOptions(TH1 *h, double max, std::string title, int ndiv, Color_t col, bool fill)
Set common options for a TLegend.
T max(const caf::Proxy< T > &a, T b)
void split(double tt, double *fr)
TString fin
Definition: Style.C:24
const XML_Char * name
Definition: expat.h:151
int nBins
Definition: plotROC.py:16
TSpline3 lo("lo", xlo, ylo, 12,"0")
Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &syst) const override
std::string fDet
Definition: PPFXHelper.h:104
TLatex * DrawLatexLines(double x, double y, std::string sPOT)
Draw common TLatex lines.
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
std::string fComponent
Definition: PPFXHelper.h:103
static std::unique_ptr< PredictionCombinePeriods > LoadFrom(TDirectory *dir, const std::string &name)
const Color_t kCosmicBackgroundColor
Definition: Style.h:41
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:209
void PlotSysts(PredictionCombinePeriods *pred, const std::vector< const ISyst * > &systs, osc::IOscCalc *calc, TDirectory *out, FILE *text, strings strs, Spectrum &sCos, const Spectrum &sOOT, Spectrum *sData)
#define MAXSPEC
Definition: NusLoadProd3.h:18
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
void OverridePOT(double newpot)
DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN&#39;T A BETTER WAY!
Definition: Spectrum.h:237
T sqrt(T number)
Definition: d0nt_math.hpp:156
virtual Spectrum PredictComponentSyst(osc::IOscCalc *calc, const SystShifts &syst, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const
Adapt the PMNS_Sterile calculator to standard interface.
string filename
Definition: shutoffs.py:106
TLegend * leg1
Definition: plot_hist.C:105
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
const Color_t kTotalMCErrorBandColor
Definition: Style.h:17
osc::OscCalcDumb calc
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
const Color_t kNumuBackgroundColor
Definition: Style.h:30
float abs(float number)
Definition: d0nt_math.hpp:39
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
fclose(fg1)
Double_t scale
Definition: plot.C:25
Charged-current interactions.
Definition: IPrediction.h:39
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
TSpline3 hi("hi", xhi, yhi, 18,"0")
std::vector< const ISyst * > getAllNusSysts()
Get a vector of all the nus group systs.
Definition: NusSysts.cxx:202
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:607
const Color_t kBeamNueBackgroundColor
Definition: Style.h:24
std::string fSystL
Definition: PPFXHelper.h:109
std::string fPOT
Definition: PPFXHelper.h:105
Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
double POT() const
Definition: Spectrum.h:231
void Preliminary()
double sigma(TH1F *hist, double percentile)
OStream cout
Definition: OStream.cxx:6
std::string fXLabel
Definition: PPFXHelper.h:110
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void SetLegendOptions(TLegend *leg)
Set common options for a TLegend.
const IPrediction * GetPrediction(int i) const
TH1 * GetHist(strings strs, int sigma)
std::string fSystType
Definition: PPFXHelper.h:107
Neutral-current interactions.
Definition: IPrediction.h:40
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
std::string fSystS
Definition: PPFXHelper.h:108
Sum MC predictions from different periods scaled according to data POT targets.
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
const Color_t kTotalMCColor
Definition: Style.h:16
T min(const caf::Proxy< T > &a, T b)
All neutrinos, any flavor.
Definition: IPrediction.h:26
const double kSecondAnaPOT
Definition: Exposures.h:87
Float_t w
Definition: plot.C:20
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:234
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:78
A helper structure to contain a group of string for plotting.
Definition: PPFXHelper.h:101
void FDDataMCSystBandAna()
std::string fSample
Definition: PPFXHelper.h:106
Spectrum Predict(osc::IOscCalc *calc) const override