Plots.cxx
Go to the documentation of this file.
2 
5 #include "CAFAna/Core/Ratio.h"
6 #include "CAFAna/Core/Spectrum.h"
8 
9 #include "Utilities/func/MathUtil.h"
10 
11 #include "TCanvas.h"
12 #include "TFeldmanCousins.h"
13 #include "TGraph.h"
14 #include "TGraphAsymmErrors.h"
15 #include "TH1.h"
16 #include "TH2.h"
17 #include "THStack.h"
18 #include "THLimitsFinder.h"
19 #include "TLegend.h"
20 #include "TLine.h"
21 #include "TList.h"
22 #include "TMarker.h"
23 #include "TPaveText.h"
24 #include "TText.h"
25 
26 #include <algorithm>
27 #include <iostream>
28 #include <cmath>
29 #include <ctime>
30 #include <functional>
31 #include <memory>
32 
33 namespace ana
34 {
35  //----------------------------------------------------------------------
36  TH1* DataMCComparison(const Spectrum& data, const Spectrum& mc, EBinType bintype)
37  {
38  TH1* ret = 0;
39 
40  TH1* hMC = mc.ToTH1(data.POT(), kTotalMCColor, kSolid, kPOT, bintype);
41 
42  TH1* hData = data.ToTH1(data.POT(), kBlack, kSolid, kPOT, bintype);
43  hData->Sumw2();
44  hData->SetMarkerStyle(kFullCircle);
45 
46  // Need to draw the biggest thing first to get correct axis limits
47  if(hMC->GetMaximum() > hData->GetMaximum()+sqrt(hData->GetMaximum())){
48  ret = hMC;
49  hMC->Draw("hist");
50  hData->Draw("ep same");
51  }
52  else{
53  ret = hData;
54  hData->Draw("ep");
55  hMC->Draw("hist same");
56  hData->Draw("ep same"); // Data always goes on top
57  }
58 
59  gPad->Update();
60 
61  return ret;
62  }
63 
64  //----------------------------------------------------------------------
66  const IPrediction* mc,
68  const SystShifts & shifts,
69  EBinType bintype)
70  {
71  return DataMCComparison(data, shifts.IsNominal() ? mc->Predict(calc) : mc->PredictSyst(calc, shifts), bintype);
72  }
73 
74  //----------------------------------------------------------------------
76  {
77  TH1* ret = 0;
78 
79  TH1* hMC = mc.ToTH1(mc.POT());
80  hMC->SetLineColor(kTotalMCColor);
81 
82  TH1* hData = data.ToTH1(data.POT());
83  hData->Sumw2();
84  hData->SetMarkerStyle(kFullCircle);
85 
86  hMC->Scale(hData->Integral()/hMC->Integral());
87 
88  // Need to draw the biggest thing first to get correct axis limits
89  if(hMC->GetMaximum() > hData->GetMaximum()+sqrt(hData->GetMaximum())){
90  ret = hMC;
91  hMC->Draw("hist");
92  hData->Draw("ep same");
93  }
94  else{
95  ret = hData;
96  hData->Draw("ep");
97  hMC->Draw("hist same");
98  hData->Draw("ep same"); // Data always goes on top
99  }
100 
101  gPad->Update();
102 
103  return ret;
104  }
105 
106  //----------------------------------------------------------------------
108  const IPrediction* mc,
110  {
111  return DataMCComparisonAreaNormalized(data, mc->Predict(calc));
112  }
113 
114  //----------------------------------------------------------------------
116  const IPrediction* mc,
118  {
119  TH1* ret = 0;
120 
121  TH1* hMC = mc->Predict(calc).ToTH1(data.POT());
122  hMC->SetLineColor(kTotalMCColor);
123 
124  TH1* hNC = mc->PredictComponent(calc,
126  Current::kNC,
127  Sign::kBoth).ToTH1(data.POT());
128  hNC->SetLineColor(kNCBackgroundColor);
129 
130  TH1* hCC = mc->PredictComponent(calc,
132  Current::kCC,
133  Sign::kBoth).ToTH1(data.POT());
134  hCC->SetLineColor(kNumuBackgroundColor);
135 
136  TH1* hBeam = mc->PredictComponent(calc,
138  Current::kCC,
139  Sign::kBoth).ToTH1(data.POT());
140  hBeam->SetLineColor(kBeamNueBackgroundColor);
141 
142  TH1* hData = data.ToTH1(data.POT());
143  hData->SetMarkerStyle(kFullCircle);
144 
145  // Need to draw the biggest thing first to get correct axis limits
146  if(hMC->GetMaximum() > hData->GetMaximum()+sqrt(hData->GetMaximum())){
147  ret = hMC;
148  hMC->Draw("hist");
149  hData->Draw("ep same");
150  }
151  else{
152  ret = hData;
153  hData->Draw("ep");
154  hMC->Draw("hist same");
155  hData->Draw("ep same"); // Data always goes on top
156  }
157 
158  hNC->Draw("hist same");
159  hCC->Draw("hist same");
160  hBeam->Draw("hist same");
161 
162  gPad->Update();
163 
164  return ret;
165  }
166 
167  //----------------------------------------------------------------------
168  void DataMCRatio(const Spectrum& data,
169  const IPrediction* mc,
171  double miny, double maxy)
172  {
173  DataMCRatio(data, mc->Predict(calc), miny, maxy);
174  }
175 
176  //----------------------------------------------------------------------
177  void DataMCRatio(const Spectrum& data,
178  const Spectrum& mc,
179  double miny, double maxy)
180  {
181  Ratio ratio(data, mc);
182  TH1* h = ratio.ToTH1();
183  h->GetYaxis()->SetTitle("Data / MC");
184  h->GetYaxis()->SetRangeUser(miny, maxy);
185  h->Draw();
186 
187  TLine* one = new TLine(h->GetXaxis()->GetXmin(), 1,
188  h->GetXaxis()->GetXmax(), 1);
189  one->SetLineWidth(2);
190  one->SetLineStyle(7);
191  one->Draw();
192 
193  gPad->Update();
194  }
195 
196  //----------------------------------------------------------------------
198  const IPrediction* mc,
200  double miny, double maxy)
201  {
202  DataMCAreaNormalizedRatio(data, mc->Predict(calc), miny, maxy);
203  }
204 
205  //----------------------------------------------------------------------
207  const Spectrum& mc,
208  double miny, double maxy)
209  {
210  Spectrum mcScaled = mc;
211  mcScaled.OverridePOT(mcScaled.POT() * mc.Integral(1) / data.Integral(1));
212 
213  DataMCRatio(data, mcScaled, miny, maxy);
214  }
215 
216  //----------------------------------------------------------------------
217  void RatioPlot(const Spectrum& data,
218  const Spectrum& expected,
219  const Spectrum& fit,
220  double miny, double maxy)
221  {
222  Ratio fitRatio(fit, expected);
223  Ratio dataRatio(data, expected);
224 
225  TH1* h = fitRatio.ToTH1();
226  h->GetYaxis()->SetTitle("Ratio to expectation");
227  h->GetYaxis()->SetRangeUser(miny, maxy);
228  h->SetLineColor(kTotalMCColor);
229  h->Draw("][");
230 
231  h = dataRatio.ToTH1();
232  h->SetMarkerStyle(kFullCircle);
233  h->Draw("ep same");
234 
235  TLine* one = new TLine(h->GetXaxis()->GetXmin(), 1,
236  h->GetXaxis()->GetXmax(), 1);
237  one->SetLineWidth(2);
238  one->SetLineStyle(7);
239  one->Draw();
240 
241  gPad->Update();
242  }
243 
244  //----------------------------------------------------------------------
246  IPrediction* pred,
247  std::vector<const ISyst*> systs,
248  std::vector<TH1*> &hUps,
249  std::vector<TH1*> &hDns,
250  double pot)
251  {
252  for (const ISyst* syst: systs) {
253  SystShifts shifts;
254 
255  shifts.SetShift(syst, +1);
256  TH1* hUp = pred->PredictSyst(calc, shifts).ToTH1(pot);
257 
258  shifts.SetShift(syst, -1);
259  TH1* hDn = pred->PredictSyst(calc, shifts).ToTH1(pot);
260 
261  hUps.push_back(hUp);
262  hDns.push_back(hDn);
263  }
264 
265  return;
266  }
267 
268  //----------------------------------------------------------------------
270  IPrediction* pred,
271  std::vector<const ISyst*> systs,
272  const SystShifts &bfshifts,
273  std::vector<TH1*> &hUps,
274  std::vector<TH1*> &hDns,
275  double pot)
276  {
277  for (const ISyst* syst: systs) {
278  SystShifts shifts;
279 
280  shifts.SetShift(syst, +1);
281  TH1* hUp = pred->PredictSyst(calc, shifts).ToTH1(pot);
282  hUp->Add( pred->PredictSyst(calc, SystShifts::Nominal()).ToTH1(pot), -1 ); // remove the nominal pred
283  hUp->Add( pred->PredictSyst(calc, bfshifts).ToTH1(pot), +1 ); // add on the BF shifted pred
284 
285  shifts.SetShift(syst, -1);
286  TH1* hDn = pred->PredictSyst(calc, shifts).ToTH1(pot);
287  hDn->Add( pred->PredictSyst(calc, SystShifts::Nominal()).ToTH1(pot), -1 ); // remove the nominal pred
288  hDn->Add( pred->PredictSyst(calc, bfshifts).ToTH1(pot), +1 ); // add on the BF shifted pred
289 
290  hUps.push_back(hUp);
291  hDns.push_back(hDn);
292  }
293 
294  return;
295  }
296 
297 
298  //----------------------------------------------------------------------
300  const std::vector<const ISyst*>& systs,
302  double pot,
303  int col, int errCol, float headroom,
304  bool newaxis,
305  EBinType bintype)
306  {
307 
308  Spectrum nom = pred->Predict(calc);
309 
310  std::vector<Spectrum> ups, dns;
311 
312  for(const ISyst* syst: systs){
313  SystShifts shifts;
314  shifts.SetShift(syst, +1);
315  ups.push_back(pred->PredictSyst(calc, shifts));
316  shifts.SetShift(syst, -1);
317  dns.push_back(pred->PredictSyst(calc, shifts));
318  }
319 
320  return PlotWithSystErrorBand(nom, ups, dns, pot, col, errCol, headroom, newaxis, bintype);
321 
322 
323  }
324 
325  //----------------------------------------------------------------------
326 
327  TGraphAsymmErrors* PlotWithSystErrorBand(const Spectrum& nominal,
328  const std::vector<Spectrum>& upShifts,
329  const std::vector<Spectrum>& downShifts,
330  double pot, int col, int errCol,
331  float headroom, bool newaxis,
332  EBinType bintype)
333  {
334 
335  TH1* nom = nominal.ToTH1(pot, kBlack, kSolid, kPOT, bintype);
336 
337  std::vector<TH1*> ups, dns;
338 
339  for(const auto& upShift:upShifts) ups.push_back(upShift.ToTH1(pot, kBlack, kSolid, kPOT, bintype));
340  for(const auto& downShift:downShifts) dns.push_back(downShift.ToTH1(pot, kBlack, kSolid, kPOT, bintype));
341 
342  return PlotWithSystErrorBand(nom, ups, dns, col, errCol, headroom, newaxis);
343  }
344 
345  //----------------------------------------------------------------------
346 
347  TGraphAsymmErrors* PlotWithSystErrorBand(TH1*& nom,
348  std::vector<TH1*>& ups,
349  std::vector<TH1*>& dns,
350  int col, int errCol,
351  float headroom, bool newaxis)
352  {
353  if(col == -1){
354  col = kTotalMCColor;
355  errCol = kTotalMCErrorBandColor;
356  }
357  else if(errCol == -1) errCol = col-7; // hopefully a lighter version
358 
359  nom->SetLineColor(col);
360  nom->GetXaxis()->CenterTitle();
361  nom->GetYaxis()->CenterTitle();
362  if(newaxis) nom->Draw("hist ]["); // Set the axes up
363 
364  double yMax = nom->GetBinContent(nom->GetMaximumBin());
365 
366  TGraphAsymmErrors* g = new TGraphAsymmErrors;
367 
368  for(int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
369  const double y = nom->GetBinContent(binIdx);
370  g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y);
371 
372  const double w = nom->GetXaxis()->GetBinWidth(binIdx);
373 
374  double errUp = 0;
375  double errDn = 0;
376 
377  for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
378  double hi = ups[systIdx]->GetBinContent(binIdx)-y;
379  double lo = dns[systIdx]->GetBinContent(binIdx)-y;
380 
381  // If both shifts are in the same direction use the larger
382  double min = std::min(hi,lo);
383  double max = std::max(hi,lo);
384  if(max < 0) max=0;
385  if(min > 0) min=0;
386 
387  errUp += max*max;
388  errDn += min*min;
389  } // end for systIdx
390 
391  g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
392  } // end for i
393 
394  g->SetFillColor(errCol);
395  g->Draw("e2 same");
396  g->GetYaxis()->SetRangeUser(0.00001, headroom*yMax);
397  nom->GetYaxis()->SetRangeUser(0.0001, headroom*yMax);
398 
399  nom->Draw("hist ][ same");
400 
401  for(TH1* up: ups) delete up;
402  for(TH1* dn: dns) delete dn;
403  return g;
404  }
405 
406  //----------------------------------------------------------------------
408  const std::vector<const ISyst*>& systs,
410  double pot,
411  int col, int errCol, float headroom,
412  bool newaxis, EBinType bintype)
413  {
414 
415  Spectrum nom = pred->Predict(calc);
416 
417  std::vector<Spectrum> ups, dns;
418 
419  for(const ISyst* syst: systs){
420  SystShifts shifts;
421  shifts.SetShift(syst, +1);
422  ups.push_back(pred->PredictSyst(calc, shifts));
423  shifts.SetShift(syst, -1);
424  dns.push_back(pred->PredictSyst(calc, shifts));
425  }
426 
427  PlotWithAreaSystErrorBand(nom, ups, dns, pot, col, errCol, headroom, newaxis, bintype);
428 
429 
430  }
431 
432  //----------------------------------------------------------------------
433  void PlotWithAreaSystErrorBand(const Spectrum& nominal,
434  std::vector<Spectrum> upShifts,
435  std::vector<Spectrum> downShifts,
436  double pot, int col, int errCol,
437  float headroom, bool newaxis,
438  EBinType bintype)
439  {
440  double nomIntegral = nominal.Integral(pot);
441  if (nomIntegral > 0)
442  {
443  // ugh so much boilerplate
444  for (auto & collection : std::vector<std::reference_wrapper<std::vector<Spectrum>>>{ upShifts, downShifts })
445  {
446  for (Spectrum & shiftSpec : collection.get())
447  {
448  // this is a dirty trick, but these Spectra are going to be thrown out once this method exits anyway
449  shiftSpec.OverridePOT(shiftSpec.POT() * shiftSpec.Integral(pot) / nomIntegral);
450  }
451  }
452  }
453 
454  PlotWithSystErrorBand(nominal, upShifts, downShifts, pot, col, errCol, headroom, newaxis, bintype);
455  }
456 
457  //----------------------------------------------------------------------
459  const std::vector<Spectrum>& upShifts,
460  const std::vector<Spectrum>& downShifts,
461  const Spectrum& alternative,
462  const std::vector<Spectrum>& upShiftsAlt,
463  const std::vector<Spectrum>& downShiftsAlt,
464  double pot, int col1, int col2,
465  float headroom, bool newaxis,
466  EBinType bintype)
467  {
468  TH1* nom = nominal.ToTH1(pot, kBlack, kSolid, kPOT, bintype);
469  TH1* alt = alternative.ToTH1(pot, kBlack, kSolid, kPOT, bintype);
470  std::vector<TH1*> ups, dns;
471  std::vector<TH1*> upsA, dnsA;
472  for(const auto& upShift:upShifts) ups.push_back(upShift.ToTH1(pot, kBlack, kSolid, kPOT, bintype));
473  for(const auto& downShift:downShifts) dns.push_back(downShift.ToTH1(pot, kBlack, kSolid, kPOT, bintype));
474  for(const auto& upShiftAlt:upShiftsAlt) upsA.push_back(upShiftAlt.ToTH1(pot, kBlack, kSolid, kPOT, bintype));
475  for(const auto& downShiftAlt:downShiftsAlt) dnsA.push_back(downShiftAlt.ToTH1(pot, kBlack, kSolid, kPOT, bintype));
476  nom->SetLineColor(col1);
477  nom->GetXaxis()->CenterTitle();
478  nom->GetYaxis()->CenterTitle();
479  if(newaxis) nom->Draw("hist"); // Set the axes up
480  alt->SetLineColor(col2);
481  double yMax = nom->GetBinContent(nom->GetMaximumBin());
482  TGraphAsymmErrors* g = new TGraphAsymmErrors;
483  TGraphAsymmErrors* g2 = new TGraphAsymmErrors;
484  for(int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
485  const double y = nom->GetBinContent(binIdx);
486  const double y2 = alt->GetBinContent(binIdx);
487  g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y);
488  g2->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y2);
489  const double w = nom->GetXaxis()->GetBinWidth(binIdx);
490  double errUp = 0, errUp2 = 0;
491  double errDn = 0, errDn2 = 0;
492  for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
493  double hi = ups[systIdx]->GetBinContent(binIdx)-y;
494  double hi2 = upsA[systIdx]->GetBinContent(binIdx)-y2;
495  double lo = y-dns[systIdx]->GetBinContent(binIdx);
496  double lo2 = y2-dnsA[systIdx]->GetBinContent(binIdx);
497  if(lo <= 0 && hi <= 0) std::swap(lo, hi);
498  if(lo2 <= 0 && hi2 <= 0) std::swap(lo2, hi2);
499  errUp += hi*hi;
500  errDn += lo*lo;
501  errUp2 += hi2*hi2;
502  errDn2 += lo2*lo2;
503  // TODO: what happens if they're both high or both low?
504  } // end for systIdx
505  g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
506  g2->SetPointError(binIdx, w/2, w/2, sqrt(errDn2), sqrt(errUp2));
507  } // end for i
508  g2->SetFillColorAlpha(col2, 0.3);
509  g2->Draw("2");
510  // g2->GetYaxis()->SetRangeUser(0, headroom*yMax);
511  g2->GetYaxis()->SetRangeUser(0.0, 19.0);
512  g->SetFillColorAlpha(col1, 0.3);
513  g->Draw("2");
514  // g->GetYaxis()->SetRangeUser(0, headroom*yMax);
515  // nom->GetYaxis()->SetRangeUser(0, headroom*yMax);
516  g->GetYaxis()->SetRangeUser(0.0, yMax*headroom);
517  nom->GetYaxis()->SetRangeUser(0.0, yMax*headroom);
518  alt->Draw("hist same");
519  nom->Draw("hist same");
520  for(TH1* up: ups) delete up;
521  for(TH1* dn: dns) delete dn;
522  for(TH1* up2: upsA) delete up2;
523  for(TH1* dn2: dnsA) delete dn2;
524  }
525 
526 
527 
528 
529 
530  //----------------------------------------------------------------------
531  TGraphAsymmErrors* PlotWithSystErrorBand_Quant(const int quant,
532  IPrediction* pred,
533  const std::vector<const ISyst*>& systs,
535  double pot,
536  int col, int errCol,
537  float maxy,
538  bool newaxis)
539  {
540 
541  Spectrum nom = pred->Predict(calc);
542 
543  std::vector<Spectrum> ups, dns;
544 
545  for(const ISyst* syst: systs){
546  SystShifts shifts;
547  shifts.SetShift(syst, +1);
548  ups.push_back(pred->PredictSyst(calc, shifts));
549  shifts.SetShift(syst, -1);
550  dns.push_back(pred->PredictSyst(calc, shifts));
551  }
552 
553  return PlotWithSystErrorBand_Quant(quant, nom, ups, dns, pot, col, errCol, maxy, newaxis);
554 
555 
556  }
557 
558  //----------------------------------------------------------------------
559  TGraphAsymmErrors* PlotWithSystErrorBand_Quant(const int quant,
560  const Spectrum& nominal,
561  const std::vector<Spectrum>& upShifts,
562  const std::vector<Spectrum>& downShifts,
563  double pot, int col, int errCol,
564  float maxy, bool newaxis)
565  {
566 
567  TH1* nom = nominal.ToTH1(pot);
568  nom->GetYaxis()->SetTitle("Events/0.1 GeV");
569  std::vector<TH1*> ups, dns;
570 
571  for(const auto& upShift:upShifts) ups.push_back(upShift.ToTH1(pot));
572  for(const auto& downShift:downShifts) dns.push_back(downShift.ToTH1(pot));
573 
574  TGraphAsymmErrors* g = PlotWithSystErrorBand_Quant(quant, nom, ups, dns, col, errCol, maxy, newaxis);
575 
576  for(TH1* up: ups) delete up;
577  for(TH1* dn: dns) delete dn;
578 
579  return g;
580  }
581 
582  //----------------------------------------------------------------------
583  TGraphAsymmErrors* PlotWithSystErrorBand_Quant(const int quant,
584  TH1*& nom,
585  std::vector<TH1*>& ups,
586  std::vector<TH1*>& dns,
587  int col, int errCol,
588  float maxy, bool newaxis)
589  {
590  if(col == -1){
591  col = kTotalMCColor;
592  errCol = kTotalMCErrorBandColor;
593  }
594  else if(errCol == -1) errCol = col-7; // hopefully a lighter version
595 
596  nom->GetYaxis()->SetTitle("Events/0.1 GeV");
597  nom->SetLineColor(col);
598  nom->GetXaxis()->CenterTitle();
599  nom->GetYaxis()->CenterTitle();
600  if(newaxis) nom->Draw("hist ]["); // Set the axes up
601 
602  TGraphAsymmErrors* g = new TGraphAsymmErrors;
603 
604  for(int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
605  const double y = nom->GetBinContent(binIdx);
606  g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y);
607 
608  const double w = nom->GetXaxis()->GetBinWidth(binIdx);
609 
610  double errUp = 0;
611  double errDn = 0;
612 
613  for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
614  double hi = ups[systIdx]->GetBinContent(binIdx)-y;
615  double lo = dns[systIdx]->GetBinContent(binIdx)-y;
616 
617  // If both shifts are in the same direction use the larger
618  double min = std::min(hi,lo);
619  double max = std::max(hi,lo);
620  if(max < 0) max=0;
621  if(min > 0) min=0;
622 
623  errUp += max*max;
624  errDn += min*min;
625  } // end for systIdx
626 
627  g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
628  } // end for i
629 
630 
631  // setting 2x2 = 4 pads canvas
632  if(quant==1 || quant==2){
633  nom->GetXaxis()->SetLabelSize(0.);
634  }
635  if(quant==2 || quant==4){
636  nom->GetYaxis()->SetLabelSize(0.);
637  }
638  float minx = 0.99 * nom->GetXaxis()->GetXmin();
639  float maxx = 0.99 * nom->GetXaxis()->GetXmax();
640  nom->GetXaxis()->SetLabelOffset(0.01);
641  nom->GetYaxis()->SetLabelOffset(0.01);
642  nom->GetYaxis()->SetTitleOffset(0.50);
643  nom->GetXaxis()->SetTitleSize(0.0);
644  nom->GetYaxis()->SetTitleSize(0.0);
645  nom->GetXaxis()->SetTitleSize(0.0);
646  nom->GetYaxis()->SetTitleSize(0.0);
647  nom->GetYaxis()->SetRangeUser(0.00001, maxy);
648  nom->GetXaxis()->SetRangeUser(minx, maxx);
649  g->GetYaxis() ->SetRangeUser(0.00001, maxy);
650  g->GetXaxis() ->SetRangeUser(minx, maxx);
651  g->SetFillColor(errCol);
652  g->Draw("e2 same");
653  nom->Draw("hist ][ same");
654 
655  return g;
656  }
657 
658 
659 
660 
661 
662 
663  //----------------------------------------------------------------------
664  THStack* ToTHStack(const std::vector<std::pair<const Spectrum&, Color_t>>& s,
665  double pot)
666  {
667  THStack* ret = new THStack;
668  for(auto it: s){
669  TH1* h = it.first.ToTH1(pot);
670  h->SetFillColor(it.second);
671  ret->Add(h, "hist");
672  }
673  return ret;
674  }
675 
676  //----------------------------------------------------------------------
677  /// Helper for \ref AutoPlaceLegend
678  double PointDistanceToBox(double x, double y,
679  double x0, double y0, double x1, double y1)
680  {
681  // Inside
682  if(x > x0 && x < x1 && y > y0 && y < y1) return 0;
683 
684  // Corners
685  double d = util::sqr(x-x0)+util::sqr(y-y0);
686  d = std::min(d, util::sqr(x-x1)+util::sqr(y-y0));
687  d = std::min(d, util::sqr(x-x1)+util::sqr(y-y1));
688  d = std::min(d, util::sqr(x-x0)+util::sqr(y-y1));
689 
690  // Top and bottom edges
691  if(x > x0 && x < x1){
692  d = std::min(d, util::sqr(y-y0));
693  d = std::min(d, util::sqr(y-y1));
694  }
695  // Left and right
696  if(y > y0 && y < y1){
697  d = std::min(d, util::sqr(x-x0));
698  d = std::min(d, util::sqr(x-x1));
699  }
700 
701  return d;
702  }
703 
704  //----------------------------------------------------------------------
705  TLegend* AutoPlaceLegend(double dx, double dy, double yPin)
706  {
707  gPad->Update();
708 
709  // Convert requested width and height into physics coordinates
710  dx *= (gPad->GetX2()-gPad->GetX1());
711  dy *= (gPad->GetY2()-gPad->GetY1());
712 
713  // Range of axes in physics units
714  const double x0 = gPad->GetUxmin();
715  const double x1 = gPad->GetUxmax();
716  const double y0 = gPad->GetUymin();
717  const double y1 = gPad->GetUymax();
718 
719  const double X = x1-x0;
720  const double Y = y1-y0;
721 
722  double bestd = 0;
723  double bestx = 0;
724  double besty = 0;
725  // If we want to pin Y pos, set it to that now.
726  if(yPin >= 0) besty = yPin * (gPad->GetY2()-gPad->GetY1());
727 
728  for(bool fallback: {false, true}){
729  for(double x = x0+dx/2; x < x1-dx/2; x += X/50){
730  for(double y = y0+dy/2; y < y1-dy/2; y += Y/50){
731  double d = 999999;
732 
733  // Repel from edges
734  d = std::min(d, util::sqr((x-dx/2-x0)/X));
735  d = std::min(d, util::sqr((x+dx/2-x1)/X));
736  d = std::min(d, util::sqr((y-dy/2-y0)/Y));
737  d = std::min(d, util::sqr((y+dy/2-y1)/Y));
738  if(d < bestd) continue;
739 
740  TIter next(gPad->GetListOfPrimitives());
741  while(TObject* obj = next()){
742 
743  if(!obj->InheritsFrom(TH1::Class())) continue;
744  if( obj->InheritsFrom(TH2::Class())) continue;
745 
746  TH1* h = (TH1*)obj;
747 
748  for(int n = 1; n <= h->GetNbinsX(); ++n){
749  const double px = h->GetBinCenter(n);
750  const double py = h->GetBinContent(n);
751 
752  if(fallback){
753  d = std::min(d, util::sqr((px-x)/X)+util::sqr((py-y)/Y));
754  }
755  else{
756  d = std::min(d, PointDistanceToBox(px/X, py/Y,
757  (x-dx/2)/X, (y-dy/2)/Y,
758  (x+dx/2)/X, (y+dy/2)/Y));
759  }
760  if(d < bestd) break;
761  }
762  if(d < bestd) break;
763  } // end while
764 
765  if(d > bestd){
766  bestd = d;
767  bestx = x;
768  // Update Y if we're not pinning it.
769  if (yPin < 0) besty = y;
770  }
771  } // end for y
772  } // end for x
773 
774  if(bestd != 0) break; // If we always collide, have to do fallback
775  } // end for fallback
776 
777  // Convert to pad coordinates
778  const double nx = (bestx-gPad->GetX1())/(gPad->GetX2()-gPad->GetX1());
779  const double ny = (besty-gPad->GetY1())/(gPad->GetY2()-gPad->GetY1());
780 
781  const double ndx = dx/(gPad->GetX2()-gPad->GetX1());
782  const double ndy = dy/(gPad->GetY2()-gPad->GetY1());
783 
784  return new TLegend(nx-ndx/2, ny-ndy/2, nx+ndx/2, ny+ndy/2);
785  }
786 
787  //----------------------------------------------------------------------
788  void ErrorBarChart(const std::map<std::string, std::pair<double, double>>& systErrors,
789  const std::pair<double, double> & statErr,
790  const std::string & label )
791  {
792  bool doStat = statErr.first != 0 && statErr.second != 0;
793  int nbins = systErrors.size() + 1; // extra for the "total systematic" bar
794 
795  // Systematics are on top, total systematic and (optional) statistical
796  // error bar are at the bottom.
797  int firstSystBin = 1;
798  if(doStat){
799  ++nbins;
800  ++firstSystBin;
801  }
802 
803  // Sum all the systematics in qaudrature
804  std::pair<double, double> systTot = std::make_pair(0., 0.);
805  for(auto it_systPairs: systErrors)
806  {
807  systTot.first += util::sqr(it_systPairs.second.first);
808  systTot.second += util::sqr(it_systPairs.second.second);
809  }
810  systTot.first = sqrt(systTot.first);
811  systTot.second = sqrt(systTot.second);
812 
813  // x-axis range is a multiple of 5% with at least 2% padding
814  double xrange = 0;
815  while(xrange < std::max(std::max(systTot.first, systTot.second), std::max(statErr.first, statErr.second))+2) xrange += 5;
816 
817  TH2* axes = new TH2F("", (";" + label).c_str(),
818  10, -xrange, +xrange, nbins, 0, nbins);
819  axes->GetXaxis()->CenterTitle();
820 
821  // Put the systematics in a vector so we can sort them by size
822  std::vector<std::pair<std::pair<double, double>, std::string>> systList;
823  for(auto it: systErrors) systList.push_back(std::make_pair(it.second, it.first));
824  std::sort(systList.begin(),
825  systList.end(),
826  [](const auto & errNamePair1, const auto & errNamePair2)
827  {
828  return ( std::max(errNamePair1.first.first, errNamePair1.first.second)
829  < std::max(errNamePair2.first.first, errNamePair2.first.second) );
830  }
831  );
832 
833  // Label the y axis
834  TAxis* yax = axes->GetYaxis();
835  for(unsigned int i = 0; i < systList.size(); ++i)
836  yax->SetBinLabel(i+firstSystBin+1, systList[i].second.c_str());
837 
838  yax->SetBinLabel(firstSystBin, "Total syst. error");
839  if(doStat) yax->SetBinLabel(1, "Statistical error");
840 
841  // Make the labels approximately match the x axis title
842  yax->SetLabelSize(.06);
843 
844  axes->Draw();
845 
846  // Blue boxes for each systematic
847  for(unsigned int i = 0; i < systList.size(); ++i){
848  TBox* box = new TBox(-systList[i].first.first, i+firstSystBin+.2,
849  +systList[i].first.second, i+firstSystBin+.8);
850  std::cout<<"ERROR = (" << systList[i].first.first << ", " << systList[i].first.second << ")" << std::endl;
851  box->SetFillColor(kAzure-8);
852  box->Draw();
853  }
854 
855  // Total systematic error
856  TBox* syst = new TBox(-systTot.first, firstSystBin-.8,
857  +systTot.second, firstSystBin-.2);
858  std::cout<<"TOTAL SYSTEMATIC ERROR = (" << systTot.first << ", " << systTot.second << ")" << std::endl;
859  syst->SetFillColor(kAzure-8);
860  syst->Draw();
861 
862  // Red box for statistical error
863  if(doStat){
864  TBox* stat = new TBox(-statErr.first, .2, +statErr.second, .8);
865  stat->SetFillColor(kRed-7);
866  stat->Draw();
867  }
868 
869  // Separate individual systematics from the total and statistical errors
870  TLine* div = new TLine(-xrange, firstSystBin, +xrange, firstSystBin);
871  div->SetLineWidth(2);
872  div->Draw();
873 
874  // Vertical line marking the symmetry point
875  TLine* zero = new TLine(0, 0, 0, nbins);
876  zero->SetLineStyle(7);
877  zero->SetLineWidth(2);
878  zero->Draw();
879 
880  // Leave enough space for the systematic labels
881  gPad->SetLeftMargin(.25);
882 
883  }
884 
885  //----------------------------------------------------------------------
886  void CountingExperimentErrorBarChart(const std::map<std::string, double>& systbars, double statErr, bool bkgdOrSig, bool shortchart)
887  {
888  std::map<std::string, double> systs = systbars; // breaks shortchart option (shortchart) ? SumNueSecondAnaSysts(systbars) : systbars;
889 
890  std::map<std::string, std::pair<double, double>> systVals;
891  for (const auto & systPair : systbars)
892  systVals[systPair.first] = std::make_pair(systPair.second, systPair.second);
893  auto statVals = std::make_pair(statErr, statErr);
894 
895  std::string label = std::string(bkgdOrSig ? "Signal" : "Background") + " uncertainty (%)";
896  ErrorBarChart(systVals, statVals, label);
897  }
898 
899  //----------------------------------------------------------------------
900  TGraphAsymmErrors* GraphWithPoissonErrors(const TH1* h, bool noErrorsXaxis, bool drawEmptyBins)
901  {
902  TGraphAsymmErrors* gr = new TGraphAsymmErrors(h);
903 
904  TFeldmanCousins fc(0.6827);//1 sigma
905 
906  for(int i = 0; i < h->GetNbinsX(); ++i){
907  double x, y;
908  gr->GetPoint(i, x, y);
909 
910  if ( drawEmptyBins || y!=0 ){
911  if ( y < 50 ) gr->SetPointEYlow(i, y-fc.CalculateLowerLimit(y,0));
912  else gr->SetPointEYlow(i,sqrt(y));
913  if ( y < 30 ) gr->SetPointEYhigh(i,fc.CalculateUpperLimit(y,0)-y);
914  else gr->SetPointEYhigh(i,sqrt(y));
915  }
916 
917  if(noErrorsXaxis){
918  gr->SetPointEXlow(i, 0);
919  gr->SetPointEXhigh(i, 0);
920  } // Do not use bin width as X-error for points
921  }
922  gr->SetMarkerStyle(kFullCircle);
923  gr->SetMarkerColor(1);
924  gr->SetMarkerSize(1);
925  gr->SetLineWidth(2);
926 
927 
928  return gr;
929  }
930 
931 
932  //----------------------------------------------------------------------
933  TGraph* graphAsymmErrorScaled(TH1* histScaled, TH1* hist, double overallScale){
934 
935  for(int bin = 1; bin <= hist->GetNbinsX(); bin++) {
936  if(hist->GetBinContent(bin)==0){
937  hist->SetBinContent(bin,0.001);
938  histScaled->SetBinContent(bin,0.001);
939  }
940  }
941 
942  TFeldmanCousins fc(0.6827); // to get correct poisson error bars
943  TGraphAsymmErrors *datapoisson = new TGraphAsymmErrors(histScaled);
944 
945  for(int bin = 1; bin <= hist->GetNbinsX(); bin++) {
946  //default numu: scale by bin 0.1/width
947  double binWidth = hist->GetBinWidth(bin);
948  double scaleFactor = overallScale / binWidth;
949  double y = hist->GetBinContent(bin);
950 
951  // since this is data it should be a whole number
952  // fc Upper/Lower limit is finicky eg if y = 3.01
953  // the returned value is the same as when y = 4
954  y = round(y);
955 
956  double errup = 0, errdn = 0;
957  if(y < 30) errup = (fc.CalculateUpperLimit(y,0) - y);
958  else errup = sqrt(y);
959  if(y < 50) errdn = (y - fc.CalculateLowerLimit(y,0));
960  else errdn = sqrt(y);
961 
962  datapoisson->SetPointEYhigh(bin-1, scaleFactor * errup);
963  datapoisson->SetPointEYlow (bin-1, scaleFactor * errdn);
964  }
965  datapoisson->SetMarkerStyle(kFullCircle);
966  datapoisson->SetMarkerColor(1);
967  datapoisson->SetMarkerSize(1);
968  datapoisson->SetLineWidth(2);
969 
970  return datapoisson;
971 
972  }
973 
974  //----------------------------------------------------------------------
975  TGraph* graphAsymmErrorWithBkgScaled(TH1* data, TH1* bkgd, double overallScale){
976 
977  for(int bin = 1; bin <= data->GetNbinsX(); bin++) {
978  if(data->GetBinContent(bin)==0){
979  data->SetBinContent(bin,0.001);
980  }
981  }
982 
983  TFeldmanCousins fc(0.6827); // to get correct poisson error bars
984  TGraphAsymmErrors *gr = new TGraphAsymmErrors(data);
985 
986  for(int bin = 1; bin <= data->GetNbinsX(); bin++) {
987  //default numu: scale by bin 0.1/width
988  double width = data->GetBinWidth(bin);
989  double scale = width / overallScale; // this is to reverse scaling
990 
991  double x, y;
992  gr->GetPoint(bin-1, x, y);
993  double bkg = bkgd->GetBinContent(bin);
994  y *= scale;
995  bkg *= scale;
996 
997  // since this is data it should be a whole number
998  // fc Upper/Lower limit is finicky eg if y = 3.01
999  // the returned value is the same as when y = 4
1000  y = round(y);
1001 
1002  double errup = 0, errdn = 0;
1003  if(y < 30) errup = (fc.CalculateUpperLimit(y,bkg) - y);
1004  else errup = sqrt(y);
1005  if(y < 50) errdn = (y - fc.CalculateLowerLimit(y,bkg));
1006  else errdn = sqrt(y);
1007 
1008  gr->SetPointEYhigh(bin-1, 1/scale * errup);
1009  gr->SetPointEYlow (bin-1, 1/scale * errdn);
1010  }
1011 
1012  gr->SetMarkerStyle(kFullCircle);
1013  gr->SetMarkerColor(1);
1014  gr->SetMarkerSize(1);
1015  gr->SetLineWidth(2);
1016 
1017  return gr;
1018 
1019  }
1020 
1021  //----------------------------------------------------------------------
1022  TGraph* RatioAsymmError(TH1* data, TH1* pred)
1023  {
1024  TGraphAsymmErrors* gRat = GraphWithPoissonErrors(data);
1025  TFeldmanCousins fc(0.6827); // to get correct poisson error bars
1026 
1027  // Hist underflow bin is 0 so index from 1
1028  // Graphs start at 0 though so will offset points
1029  for (int bin = 1; bin <= data->GetNbinsX(); ++bin) {
1030  double x, y;
1031  gRat->GetPoint(bin-1, x, y);
1032  double den = pred->GetBinContent(bin);
1033 
1034  // since this is data it should be a whole number
1035  // fc Upper/Lower limit is finicky eg if y = 3.01
1036  // the returned value is the same as when y = 4
1037  // Also negatives are weird
1038  y = round(y);
1039  if (y <= 0) y = 0.01;
1040 
1041  double errup = 0, errdn = 0;
1042  if (y < 30) errup = (fc.CalculateUpperLimit(y, 0) - y) / den;
1043  else errup = sqrt(y) / den;
1044  if (y < 50) errdn = (y - fc.CalculateLowerLimit(y, 0)) / den;
1045  else errdn = sqrt(y) / den;
1046 
1047  gRat->SetPoint(bin-1, x, y / den);
1048  gRat->SetPointEYhigh(bin-1, errup);
1049  gRat->SetPointEYlow(bin-1, errdn);
1050  }
1051  gRat->SetMarkerStyle(kFullCircle);
1052  gRat->SetMarkerColor(1);
1053  gRat->SetMarkerSize(1);
1054  gRat->SetLineWidth(2);
1055 
1056  return gRat;
1057  }
1058 
1059  //----------------------------------------------------------------------
1060  TGraph* RatioAsymmErrorScaled(TH1* data, TH1* pred, double overallScale)
1061  {
1062  TGraphAsymmErrors* gRat = GraphWithPoissonErrors(data);
1063  TFeldmanCousins fc(0.6827); // to get correct poisson error bars
1064 
1065  // Hist underflow bin is 0 so index from 1
1066  // Graphs start at 0 though so will offset points
1067  for (int bin = 1; bin <= data->GetNbinsX(); ++bin) {
1068  double width = data->GetBinWidth(bin);
1069  double scale = width / overallScale; // this is to reverse scaling
1070 
1071  double x, y;
1072  gRat->GetPoint(bin-1, x, y);
1073  double den = pred->GetBinContent(bin);
1074  y *= scale;
1075  den *= scale;
1076 
1077  // since this is data it should be a whole number
1078  // fc Upper/Lower limit is finicky eg if y = 3.01
1079  // the returned value is the same as when y = 4
1080  // Also negatives are weird
1081  y = round(y);
1082  if (y <= 0) y = 0.01;
1083 
1084  double errup = 0, errdn = 0;
1085  if (y < 30) errup = (fc.CalculateUpperLimit(y, 0) - y) / den;
1086  else errup = sqrt(y) / den;
1087  if (y < 50) errdn = (y - fc.CalculateLowerLimit(y, 0)) / den;
1088  else errdn = sqrt(y) / den;
1089 
1090  gRat->SetPoint(bin-1, x, y / den);
1091  gRat->SetPointEYhigh(bin-1, errup);
1092  gRat->SetPointEYlow (bin-1, errdn);
1093  }
1094  gRat->SetMarkerStyle(kFullCircle);
1095  gRat->SetMarkerColor(1);
1096  gRat->SetMarkerSize(1);
1097  gRat->SetLineWidth(2);
1098 
1099  return gRat;
1100  }
1101 
1102  //----------------------------------------------------------------------
1103  TGraph* RatioAsymmErrorWithBkg(TH1* data, TH1* pred, TH1* bkgd)
1104  {
1105  TGraphAsymmErrors* gRat = GraphWithPoissonErrors(data);
1106  TFeldmanCousins fc(0.6827); // to get correct poisson error bars
1107 
1108  // Hist underflow bin is 0 so index from 1
1109  // Graphs start at 0 though so will offset points
1110  for (int bin = 1; bin <= data->GetNbinsX(); ++bin) {
1111  double x, y;
1112  gRat->GetPoint(bin-1, x, y);
1113  double den = pred->GetBinContent(bin);
1114  double bkg = bkgd->GetBinContent(bin);
1115 
1116  // since this is data it should be a whole number
1117  // fc Upper/Lower limit is finicky eg if y = 3.01
1118  // the returned value is the same as when y = 4
1119  // Also negatives are weird
1120  y = round(y);
1121  if (y <= 0) y = 0.01;
1122 
1123  double errup = 0, errdn = 0;
1124  if (y < 30) errup = (fc.CalculateUpperLimit(y, bkg) - y) / den;
1125  else errup = sqrt(y) / den;
1126  if (y < 50) errdn = (y - fc.CalculateLowerLimit(y, bkg)) / den;
1127  else errdn = sqrt(y) / den;
1128 
1129  gRat->SetPoint(bin-1, x, y / den);
1130  gRat->SetPointEYhigh(bin-1, errup);
1131  gRat->SetPointEYlow (bin-1, errdn);
1132  }
1133  gRat->SetMarkerStyle(kFullCircle);
1134  gRat->SetMarkerColor(1);
1135  gRat->SetMarkerSize(1);
1136  gRat->SetLineWidth(2);
1137 
1138  return gRat;
1139  }
1140 
1141 
1142  //----------------------------------------------------------------------
1143  TGraph* RatioAsymmErrorWithBkgScaled(TH1* data, TH1* pred, TH1* bkgd, double overallScale)
1144  {
1145  TGraphAsymmErrors* gRat = GraphWithPoissonErrors(data);
1146  TFeldmanCousins fc(0.6827); // to get correct poisson error bars
1147 
1148  // Hist underflow bin is 0 so index from 1
1149  // Graphs start at 0 though so will offset points
1150  for (int bin = 1; bin <= data->GetNbinsX(); ++bin) {
1151  double width = data->GetBinWidth(bin);
1152  double scale = width / overallScale; // this is to reverse scaling
1153 
1154  double x, y;
1155  gRat->GetPoint(bin-1, x, y);
1156  double den = pred->GetBinContent(bin);
1157  double bkg = bkgd->GetBinContent(bin);
1158 
1159  y *= scale;
1160  den *= scale;
1161  bkg *= scale;
1162 
1163  // since this is data it should be a whole number
1164  // fc Upper/Lower limit is finicky eg if y = 3.01
1165  // the returned value is the same as when y = 4
1166  // Also negatives are weird
1167  y = round(y);
1168  if (y <= 0) y = 0.01;
1169 
1170  double errup = 0, errdn = 0;
1171  if (y < 30) errup = (fc.CalculateUpperLimit(y, bkg) - y) / den;
1172  else errup = sqrt(y) / den;
1173  if (y < 50) errdn = (y - fc.CalculateLowerLimit(y, bkg)) / den;
1174  else errdn = sqrt(y) / den;
1175 
1176  gRat->SetPoint(bin-1, x, y / den);
1177  gRat->SetPointEYhigh(bin-1, errup);
1178  gRat->SetPointEYlow (bin-1, errdn);
1179  }
1180  gRat->SetMarkerStyle(kFullCircle);
1181  gRat->SetMarkerColor(1);
1182  gRat->SetMarkerSize(1);
1183  gRat->SetLineWidth(2);
1184 
1185  return gRat;
1186  }
1187 
1188  //----------------------------------------------------------------------
1189  TGraph* ShadeBetweenHistograms(TH1* hmin, TH1* hmax)
1190  {
1191  int n = hmin->GetNbinsX();
1192  TGraph* gr = new TGraph(4*n);
1193 
1194  for(int i=1; i<=n; i++)
1195  {
1196  double xdown = hmax->GetBinLowEdge(i);
1197  double xup = hmax->GetBinLowEdge(i+1);
1198  double ydown = hmin->GetBinContent(i);
1199  double yup = hmax->GetBinContent(i);
1200 
1201  gr->SetPoint(2*(i-1), xdown, ydown);
1202  gr->SetPoint(2*(i-1)+1, xup, ydown);
1203  gr->SetPoint(4*n-2*i, xup, yup);
1204  gr->SetPoint(4*n-2*i+1, xdown, yup);
1205  }
1206 
1207  gr->SetFillColor(hmin->GetLineColor() - 9); // In principle this is lighter
1208  gr->SetLineColor(hmin->GetLineColor() - 9); // In case one does Draw("lf")
1209 
1210  return gr;
1211  }
1212  //----------------------------------------------------------------------
1213  TGraphAsymmErrors * ProfileQuantile(const TH2 * hist,
1214  const std::string & axisName,
1215  const std::string & graphName,
1216  const std::pair<double, double> & quantileDivisions)
1217  {
1218  if (hist->GetDimension() != 2)
1219  throw std::runtime_error(Form("Can't profile a histogram with other than 2 dimensions. Yours is a %d-D histogram...", hist->GetDimension()));
1220 
1221  TAxis const* axis = nullptr;
1222  std::function<TH1D*(const char *, Int_t, Int_t, Option_t*)> projectionMethod;
1223  if (axisName == "x" || axisName == "X")
1224  {
1225  axis = hist->GetXaxis();
1226  projectionMethod = std::bind(&TH2::ProjectionY, hist, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4);
1227  }
1228  else if (axisName == "y" || axisName == "Y")
1229  {
1230  axis = hist->GetYaxis();
1231  projectionMethod = std::bind(&TH2::ProjectionX, hist, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4);
1232  }
1233  else
1234  throw std::runtime_error( Form("Can't profile onto unknown axis: '%s'", axisName.c_str()) );
1235 
1236 
1237  std::vector<double> points_x, points_y, errors_x_low, errors_x_high, errors_y_low, errors_y_high;
1238 
1239  double quantiles[2] = {0, 0};
1240  double _quantileDivisions[2] = {quantileDivisions.first, quantileDivisions.second};
1241  for (int binNum = 1; binNum <= axis->GetNbins(); binNum++) // remember, 0 is underflow and Nbins+1 is overflow...
1242  {
1243  std::unique_ptr<TH1D> proj ( projectionMethod(Form("%s_tmp_%d", hist->GetName(), int(std::time(nullptr))), binNum, binNum, "") ); // randomish name, just to be safe
1244 
1245  double y = proj->GetMean();
1246  points_x.push_back(axis->GetBinCenter(binNum));
1247  points_y.push_back(y);
1248 
1249  // for now, make the errors for an empty projection 0
1250  unsigned int n_qs = 0;
1251  if (proj->Integral() == 0)
1252  {
1253  n_qs = 2;
1254  quantiles[0] = quantiles[1] = 0;
1255  }
1256  else
1257  n_qs = proj->GetQuantiles(2, quantiles, _quantileDivisions);
1258  if (n_qs != 2)
1259  throw std::runtime_error( Form("GetQuantiles() didn't compute all the quantiles in HistoTools.ProfileQuantile(). I requested 2 quantiles, but got %d of them...", n_qs) );
1260 
1261  double binWidth = axis->GetBinWidth(binNum);
1262  errors_x_low.push_back(binWidth/2.);
1263  errors_x_high.push_back(binWidth/2.);
1264  errors_y_low.push_back(y-quantiles[0]);
1265  errors_y_high.push_back(quantiles[1]-y);
1266  }
1267 
1268  TGraphAsymmErrors * outGraph = new TGraphAsymmErrors(
1269  points_x.size(),
1270  &points_x[0],
1271  &points_y[0],
1272  &errors_x_low[0],
1273  &errors_x_high[0],
1274  &errors_y_low[0],
1275  &errors_y_high[0]
1276  );
1277  std::string name = (graphName.size()) ? graphName : Form("%s_quantile_errors", hist->GetName());
1278  outGraph->SetName(name.c_str());
1279 
1280  return outGraph;
1281 
1282  }
1283 
1284 
1285  // Draw single BF
1286  void drawBFSingle(double bfSin, double bfDm, Color_t color, Style_t marker, double size = 1.5)
1287  {
1288 
1289  TMarker* manMarker = new TMarker(bfSin, bfDm, marker);
1290  manMarker->SetMarkerSize(size);
1291  manMarker->SetMarkerColor(color);
1292  manMarker->Draw();
1293 
1294  }
1295 
1296  // Truth for FHC; not sure about RHC
1297  void drawBFMirror(double bfSin, double bfDm, Color_t color, Style_t marker, double size = 1.5)
1298  {
1299 
1300  double mirror = bfSin;
1301  TMarker* manMarker = new TMarker(bfSin, bfDm, marker);
1302  manMarker->SetMarkerSize(size);
1303  manMarker->SetMarkerColor(color);
1304  manMarker->Draw();
1305 
1306  if(bfSin > 0.514) mirror = 0.514 - (bfSin - 0.514);
1307  if(bfSin < 0.514) mirror = 0.514 + (0.514 - bfSin );
1308  TMarker* mirrorMarker = new TMarker(mirror, bfDm, marker);
1309  mirrorMarker->SetMarkerSize(size);
1310  mirrorMarker->SetMarkerColor(color);
1311  mirrorMarker->Draw();
1312 
1313  }
1314 
1315 
1316 
1317  //----------------------------------------------------------------------
1318  void MakeHistCanvasReady_Quant(const int quant,
1319  TH1* hist,
1320  double ymax){
1321 
1322  hist->GetXaxis()->SetLabelOffset(0.01);
1323  hist->GetYaxis()->SetLabelOffset(0.01);
1324  hist->GetYaxis()->SetTitleOffset(0.50);
1325  hist->GetXaxis()->SetTitleSize(0.0);
1326  hist->GetYaxis()->SetTitleSize(0.0);
1327  hist->GetYaxis()->SetRangeUser(0.0, ymax);
1328 
1329  if(quant==1 || quant==2) hist->GetXaxis()->SetLabelSize(0.);
1330  if(quant==2 || quant==4) hist->GetYaxis()->SetLabelSize(0.);
1331  // Q3 might need some axis labels removed to prevent interference.
1332  // Rather than assuming there will be trouble we can check the
1333  // optimization that ROOT does when plottng the axes
1334  // ChangeLabel labnum -1 alters the last label
1335  if(quant == 3) {
1336  // initialize the optimized min, max, step size,
1337  // and number of steps of the labels
1338  double optmin = 0., optmax = 0., optwi = 0.;
1339  int optstp = 0;
1340 
1341  double xmin = hist->GetXaxis()->GetXmin();
1342  double xmax = hist->GetXaxis()->GetXmax();
1343  int xdiv = hist->GetXaxis()->GetNdivisions() % 100;
1344  THLimitsFinder::Optimize(xmin, xmax, xdiv, optmin, optmax, optstp, optwi);
1345 
1346  if ( (xmax - optmax) < optwi/2 )
1347  hist->GetXaxis()->ChangeLabel( -1, -1, 0);
1348 
1349  optmin = 0.; optmax = 0.; optwi = 0.; optstp = 0;
1350  double ymin = 0.;
1351  int ydiv = hist->GetYaxis()->GetNdivisions() % 100;
1352  THLimitsFinder::Optimize(ymin, ymax, ydiv, optmin, optmax, optstp, optwi);
1353 
1354  hist->GetYaxis()->SetRangeUser(ymin, ymax);
1355  if ( (ymax - optmax) < optwi/2 )
1356  hist->GetYaxis()->ChangeLabel( -1, -1, 0);
1357  }
1358  }
1359 
1360  //----------------------------------------------------------------------
1361  void PimpHist(TH1* hist, Style_t linestyle, Color_t linecolor, int linewidth, Style_t markerstyle, Color_t markercolor, double markersize){
1362  hist->SetLineColor(linecolor);
1363  hist->SetLineStyle(linestyle);
1364  hist->SetLineWidth(linewidth);
1365  hist->SetMarkerColor(markercolor);
1366  hist->SetMarkerStyle(markerstyle);
1367  hist->SetMarkerSize(markersize);
1368  }
1369 
1370 
1371  //----------------------------------------------------------------------
1372  void SplitCanvas(double ysplit, TPad*& p1, TPad*& p2)
1373  {
1374  if(!gPad) new TCanvas;
1375  p1 = new TPad("", "", 0, 0, 1, 1);
1376  p2 = new TPad("", "", 0, 0, 1, 1);
1377 
1378  p1->SetBottomMargin(ysplit);
1379  p2->SetTopMargin(1-ysplit);
1380 
1381  // Draw p1 second since it's often the more important one, that the user
1382  // would prefer to be able to interact with.
1383  for(TPad* p: {p2, p1}){
1384  p->SetFillStyle(0);
1385  p->Draw();
1386  }
1387  }
1388 
1389 
1390  //----------------------------------------------------------------------
1391  // split canvas in 4 for numu quartiles
1392  void SplitCanvasQuant(TCanvas *& canvas, TPad *& pad1, TPad *& pad2, TPad *& pad3, TPad *& pad4){
1393 
1394  canvas = new TCanvas(UniqueName().c_str(),"",960,800);
1395 
1396  pad3 = new TPad(UniqueName().c_str(),"pad3",0,0,1,1);
1397  pad3->Draw();
1398  pad3->SetTopMargin(0.5);
1399  pad3->SetBottomMargin(0.1);
1400  pad3->SetLeftMargin(0.1);
1401  pad3->SetRightMargin(0.49);
1402  pad3->SetFillStyle(0);
1403  canvas->cd();
1404 
1405  pad4 = new TPad(UniqueName().c_str(),"pad4",0,0,1,1);
1406  pad4->Draw();
1407  pad4->SetTopMargin(0.5);
1408  pad4->SetBottomMargin(0.1);
1409  pad4->SetLeftMargin(0.51);
1410  pad4->SetRightMargin(0.08);
1411  pad4->SetFillStyle(0);
1412  canvas->cd();
1413 
1414  pad1 = new TPad(UniqueName().c_str(),"pad1",0,0,1,1);// x1 y1 x2 y2
1415  pad1->Draw();
1416  pad1->SetTopMargin(0.1);
1417  pad1->SetBottomMargin(0.5);
1418  pad1->SetLeftMargin(0.1);
1419  pad1->SetRightMargin(0.49);
1420  pad1->SetFillStyle(0);
1421  canvas->cd();
1422 
1423  pad2 = new TPad(UniqueName().c_str(),"pad2",0,0,1,1);
1424  pad2->Draw();
1425  pad2->SetTopMargin(0.1);
1426  pad2->SetBottomMargin(0.5);
1427  pad2->SetLeftMargin(0.51);
1428  pad2->SetRightMargin(0.08);
1429  pad2->SetFillStyle(0);
1430  canvas->cd();
1431 
1432  }
1433 
1434 
1435  //----------------------------------------------------------------------
1436  void CenterTitles ( TH1 * histo)
1437  {
1438  histo->GetXaxis()->CenterTitle();
1439  histo->GetYaxis()->CenterTitle();
1440  histo->GetZaxis()->CenterTitle();
1441  }
1442 
1443 
1444  //----------------------------------------------------------------------
1445  bool SortSystsName(const ISyst* systA, const ISyst* systB){
1446  return ( systA->ShortName() < systB->ShortName() );
1447  }
1448 
1449 
1450  TH1* PlotSystShifts(const SystShifts & shifts, bool sortName){
1451 
1452  auto systs = shifts.ActiveSysts();
1453  int nsysts = systs.size();
1454 
1455  // sort systematics by name
1456  if(sortName) std::sort(systs.begin(), systs.end(), SortSystsName);
1457 
1458  auto h = new TH1D ("sh",";; Pull (N#sigma)",
1459  nsysts, -0.5, nsysts + 0.5);
1460 
1461  for (int systIdx = 0; systIdx < nsysts; ++systIdx){
1462  double shiftThis = shifts.GetShift(systs[systIdx]);
1463  h->SetBinContent( systIdx + 1, shiftThis);
1464 
1465  TString tempName = systs[systIdx]->ShortName();
1466  if( tempName.Contains("Neutron") ) h->GetXaxis()->SetBinLabel( systIdx + 1, "NeutronSyst2018");
1467  else h->GetXaxis()->SetBinLabel( systIdx + 1, systs[systIdx]->ShortName().c_str());
1468  }
1469 
1470  h->SetMarkerStyle(kFullCircle);
1471  CenterTitles(h);
1472  h->SetTitleOffset(3.0);
1473  h->GetXaxis()->SetNdivisions(nsysts,kFALSE);
1474  h->GetXaxis()->SetLabelSize(0.065);
1475  h->GetXaxis()->LabelsOption("v");
1476  h->GetYaxis()->SetRangeUser(-1,1);
1477  auto c1 = new TCanvas ("c1","c1",700,350);
1478  c1->SetBottomMargin(0.5);
1479  h->Draw("lp hist");
1480  return h;
1481  }
1482 
1483  //----------------------------------------------------------------------
1484  TGraph* JoinGraphs(TGraph* a, TGraph* b, int fillcol)
1485  {
1486  TGraph* ret = (TGraph*)a->Clone(UniqueName().c_str());
1487  for(int i = 0; i < b->GetN(); ++i){
1488  double x, y;
1489  b->GetPoint(i, x, y);
1490  ret->SetPoint(ret->GetN(), x, y);
1491  }
1492 
1493  ret->SetFillColor(fillcol);
1494 
1495  return ret;
1496  }
1497 
1498  //----------------------------------------------------------------------
1499  TGraph* ExtendGraphToTop(TGraph* g, int col, double xmin, double xmax, double y)
1500  {
1501  g->SetPoint(g->GetN(), xmax, y);
1502  g->SetPoint(g->GetN(), xmin, y);
1503 
1504  g->SetFillColor(col);
1505 
1506  return g;
1507  }
1508 
1509  //----------------------------------------------------------------------
1510  TPaveText* DrawBeamLabel(bool isFHC)
1511  {
1512 
1513  TPaveText *pText = new TPaveText (0.1, 0.90, 0.16, 0.94,"NDC");
1514  pText->SetFillStyle(0);
1515  pText->SetLineColor(0);
1516  pText->SetLineWidth(0);
1517  pText->SetBorderSize(1);
1518  pText->SetTextColor(kGray+2);
1519  if (isFHC) pText->AddText("#nu-beam");
1520  else pText->AddText("#bar{#nu}-beam");
1521  pText->SetTextSize(2/40.);
1522  pText->SetTextAlign(11);
1523  pText->Draw();
1524 
1525  return pText;
1526  }
1527 
1528  //----------------------------------------------------------------------
1529  TPaveText* DrawQuantLabel(int quant)
1530  {
1531  TPaveText *ret = 0;
1532  // Draw nothing for q = 0
1533  if (quant == 0) {
1534  // Draw nothing for q = 0 but it's still valid
1535  }
1536  else if (quant == 1) {
1537  ret = new TPaveText(0.37, 0.75, 0.47, 0.90, "NDC");
1538  ret->SetBorderSize(0);
1539  ret->SetFillStyle(0);
1540  TText *text = ret->AddText(TString::Format("#splitline{#bf{Quartile %d}}{best resolution}", quant));
1541  text ->SetTextAlign(22);
1542  text ->SetTextSize(0.025);
1543  ret->Draw();
1544  }
1545  else if (quant == 2) {
1546  ret = new TPaveText(0.75, 0.75, 0.85, 0.90, "NDC");
1547  ret->SetBorderSize(0);
1548  ret->SetFillStyle(0);
1549  TText *text = ret->AddText(TString::Format("#bf{Quartile %d}", quant));
1550  text ->SetTextAlign(22);
1551  text ->SetTextSize(0.025);
1552  ret->Draw();
1553  }
1554  else if (quant == 3) {
1555  ret = new TPaveText(0.37, 0.35, 0.47, 0.50, "NDC");
1556  ret->SetBorderSize(0);
1557  ret->SetFillStyle(0);
1558  TText *text = ret->AddText(TString::Format("#bf{Quartile %d}", quant));
1559  text ->SetTextAlign(22);
1560  text ->SetTextSize(0.025);
1561  ret->Draw();
1562  }
1563  else if (quant == 4) {
1564  ret = new TPaveText(0.75, 0.35, 0.85, 0.50, "NDC");
1565  ret->SetBorderSize(0);
1566  ret->SetFillStyle(0);
1567  TText *text = ret->AddText(TString::Format("#splitline{#bf{Quartile %d}}{worst resolution}", quant));
1568  text ->SetTextAlign(22);
1569  text ->SetTextSize(0.025);
1570  ret->Draw();
1571  }
1572  else {
1573  std::cout << "\n***********************************************\n"
1574  << "\n No valid quantile chosen for DrawQuantLabel \n"
1575  << "\n***********************************************\n"
1576  << std::endl;
1577  }
1578 
1579  return ret;
1580  }
1581 
1582 } // namespace
const XML_Char * name
Definition: expat.h:151
std::vector< double > quantiles(TH1D *h)
Definition: absCal.cxx:528
TSpline3 lo("lo", xlo, ylo, 12,"0")
Cuts and Vars for the 2020 FD DiF Study.
Definition: CutFlow_header.h:5
std::map< std::string, double > xmax
set< int >::iterator it
TH1 * PlotSystShifts(const SystShifts &shifts, bool sortName)
Definition: Plots.cxx:1450
std::vector< SystGroupDef > systs
Definition: syst_header.h:384
TH1D * ToTH1(Color_t col=kBlack, Style_t style=kSolid) const
Definition: Ratio.cxx:122
void MakeHistCanvasReady_Quant(const int quant, TH1 *hist, double ymax)
Definition: Plots.cxx:1318
General interface to oscillation calculators.
Definition: FwdDeclare.h:15
Float_t y1[n_points_granero]
Definition: compare.C:5
bool IsNominal() const
Definition: SystShifts.h:43
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:564
virtual Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
virtual const std::string & ShortName() const final
The name printed out to the screen.
Definition: ISyst.h:27
TH1 * DataMCComparisonComponents(const Spectrum &data, const IPrediction *mc, osc::IOscCalculator *calc)
Plot MC broken down into flavour components, overlayed with data.
Definition: Plots.cxx:115
void DataMCRatio(const Spectrum &data, const IPrediction *mc, osc::IOscCalculator *calc, double miny, double maxy)
Plot data/MC ratio for the given spectrum. Normalize MC to Data by POT.
Definition: Plots.cxx:168
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:269
double maxy
Float_t x1[n_points_granero]
Definition: compare.C:5
Float_t den
Definition: plot.C:36
(&#39;beam &#39;)
Definition: IPrediction.h:15
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
TH1 * ratio(TH1 *h1, TH1 *h2)
bool isFHC
TGraphAsymmErrors * PlotWithSystErrorBand(IPrediction *pred, const std::vector< const ISyst * > &systs, osc::IOscCalculator *calc, double pot, int col, int errCol, float headroom, bool newaxis, EBinType bintype)
Plot prediction with +/-1sigma error band.
Definition: Plots.cxx:299
TPaveText * DrawQuantLabel(int quant)
Definition: Plots.cxx:1529
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:335
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1436
TGraphAsymmErrors * GraphWithPoissonErrors(const TH1 *h, bool noErrorsXaxis, bool drawEmptyBins)
Calculate statistical errors appropriate for small Poisson numbers.
Definition: Plots.cxx:900
void SplitCanvas(double ysplit, TPad *&p1, TPad *&p2)
Split the current pad into two vertically stacked pieces, suitable for ratio plots and so on...
Definition: Plots.cxx:1372
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:753
const Color_t kTotalMCErrorBandColor
Definition: Style.h:30
static SystShifts Nominal()
Definition: SystShifts.h:34
fvar< T > round(const fvar< T > &x)
Definition: round.hpp:23
::xsd::cxx::tree::time< char, simple_type > time
Definition: Database.h:194
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
const Color_t kNumuBackgroundColor
Definition: Style.h:42
TGraph * RatioAsymmError(TH1 *data, TH1 *pred)
Definition: Plots.cxx:1022
TGraphAsymmErrors * ProfileQuantile(const TH2 *hist, const std::string &axisName, const std::string &graphName, const std::pair< double, double > &quantileDivisions)
Calculate profile with error bars corresponding to specified quantiles of a 2D distribution (by defau...
Definition: Plots.cxx:1213
virtual Spectrum Predict(osc::IOscCalculator *calc) const =0
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
virtual Spectrum PredictSyst(osc::IOscCalculator *calc, const SystShifts &syst) const
Definition: IPrediction.cxx:86
Float_t Y
Definition: plot.C:38
void PlotWithAreaSystErrorBand(IPrediction *pred, const std::vector< const ISyst * > &systs, osc::IOscCalculator *calc, double pot, int col, int errCol, float headroom, bool newaxis, EBinType bintype)
Plot prediction with +/-1sigma shape-only error band.
Definition: Plots.cxx:407
const char * label
void PimpHist(TH1 *hist, Style_t linestyle, Color_t linecolor, int linewidth, Style_t markerstyle, Color_t markercolor, double markersize)
Pimp histogram once and for all.
Definition: Plots.cxx:1361
const XML_Char const XML_Char * data
Definition: expat.h:268
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
Double_t ymax
Definition: plot.C:25
TGraphAsymmErrors * PlotWithSystErrorBand_Quant(const int quant, IPrediction *pred, const std::vector< const ISyst * > &systs, osc::IOscCalculator *calc, double pot, int col, int errCol, float maxy, bool newaxis)
Definition: Plots.cxx:531
EBinType
Definition: Utilities.h:42
void drawBFSingle(double bfSin, double bfDm, Color_t color, Style_t marker, double size=1.5)
Draw a single BF point.
Definition: Plots.cxx:1286
const XML_Char * s
Definition: expat.h:262
const int nbins
Definition: cellShifts.C:15
Double_t scale
Definition: plot.C:25
Charged-current interactions.
Definition: IPrediction.h:39
TSpline3 hi("hi", xhi, yhi, 18,"0")
double dy[NP][NC]
double dx[NP][NC]
void RatioPlot(const Spectrum &data, const Spectrum &expected, const Spectrum &fit, double miny, double maxy)
Plot data/expected, compared with fit/expected.
Definition: Plots.cxx:217
void PlotWithSystErrorBandTwoPreds(const Spectrum &nominal, const std::vector< Spectrum > &upShifts, const std::vector< Spectrum > &downShifts, const Spectrum &alternative, const std::vector< Spectrum > &upShiftsAlt, const std::vector< Spectrum > &downShiftsAlt, double pot, int col1, int col2, float headroom, bool newaxis, EBinType bintype)
Plot two different predictions with +/-1sigma shape-only error bands.
Definition: Plots.cxx:458
TH1 * DataMCComparisonAreaNormalized(const Spectrum &data, const Spectrum &mc)
Definition: Plots.cxx:75
TGraph * graphAsymmErrorScaled(TH1 *histScaled, TH1 *hist, double overallScale)
Definition: Plots.cxx:933
TPaveText * DrawBeamLabel(bool isFHC)
Put the standardized beam label in the left corner of the active canvas.
Definition: Plots.cxx:1510
Int_t col[ntarg]
Definition: Style.C:29
const double a
const Color_t kBeamNueBackgroundColor
Definition: Style.h:37
TGraph * ShadeBetweenHistograms(TH1 *hmin, TH1 *hmax)
Definition: Plots.cxx:1189
#define pot
T GetShift(const ISyst *syst) const
Float_t d
Definition: plot.C:236
TGraph * ExtendGraphToTop(TGraph *g, int col, double xmin, double xmax, double y)
Definition: Plots.cxx:1499
TH2D * histo
float bin[41]
Definition: plottest35.C:14
double POT() const
Definition: Spectrum.h:263
osc::OscCalculatorDumb calc
Represent the ratio between two spectra.
Definition: Ratio.h:8
TGraph * JoinGraphs(TGraph *a, TGraph *b, int fillcol)
Join graphs and set a fill color. Useful for contours.
Definition: Plots.cxx:1484
OStream cout
Definition: OStream.cxx:6
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
double PointDistanceToBox(double x, double y, double x0, double y0, double x1, double y1)
Helper for AutoPlaceLegend.
Definition: Plots.cxx:678
TGraph * RatioAsymmErrorWithBkg(TH1 *data, TH1 *pred, TH1 *bkgd)
Definition: Plots.cxx:1103
void CountingExperimentErrorBarChart(const std::map< std::string, double > &systbars, double statErr, bool bkgdOrSig, bool shortchart)
Make a simple plot of relative size of different errors.
Definition: Plots.cxx:886
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void GetBFSystBands(osc::IOscCalculator *calc, IPrediction *pred, std::vector< const ISyst * > systs, const SystShifts &bfshifts, std::vector< TH1 * > &hUps, std::vector< TH1 * > &hDns, double pot)
Definition: Plots.cxx:269
TGraph * graphAsymmErrorWithBkgScaled(TH1 *data, TH1 *bkgd, double overallScale)
Definition: Plots.cxx:975
void DataMCAreaNormalizedRatio(const Spectrum &data, const IPrediction *mc, osc::IOscCalculator *calc, double miny, double maxy)
Plot data/MC ratio for the given spectrum. Normalize MC to Data by area.
Definition: Plots.cxx:197
Float_t proj
Definition: plot.C:35
TGraph * RatioAsymmErrorWithBkgScaled(TH1 *data, TH1 *pred, TH1 *bkgd, double overallScale)
Definition: Plots.cxx:1143
const hit & b
Definition: hits.cxx:21
void ProjectionY(T *from, U *to)
Helper for WeightingVariable.
Neutral-current interactions.
Definition: IPrediction.h:40
void ErrorBarChart(const std::map< std::string, std::pair< double, double >> &systErrors, const std::pair< double, double > &statErr, const std::string &label)
Make a simple plot of relative size of different errors.
Definition: Plots.cxx:788
THStack * ToTHStack(const std::vector< std::pair< const Spectrum &, Color_t >> &s, double pot)
Can call like ToTHStack({{h1, kRed}, {h2, kBlue}}, pot)
Definition: Plots.cxx:664
bool SortSystsName(const ISyst *systA, const ISyst *systB)
Definition: Plots.cxx:1445
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
c1
Definition: demo5.py:24
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
Double_t ymin
Definition: plot.C:24
TPad * pad3
Definition: analysis.C:13
std::vector< const ISyst * > ActiveSysts() const
Definition: SystShifts.cxx:199
string tempName
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
Definition: Plots.cxx:705
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:29
void SplitCanvasQuant(TCanvas *&canvas, TPad *&pad1, TPad *&pad2, TPad *&pad3, TPad *&pad4)
Definition: Plots.cxx:1392
TGraph * RatioAsymmErrorScaled(TH1 *data, TH1 *pred, double overallScale)
Definition: Plots.cxx:1060
auto one()
Definition: PMNS.cxx:49
TPad * pad2
Definition: analysis.C:13
TH1 * DataMCComparison(const Spectrum &data, const Spectrum &mc, EBinType bintype)
Definition: Plots.cxx:36
All neutrinos, any flavor.
Definition: IPrediction.h:26
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
Float_t X
Definition: plot.C:38
const Color_t kNCBackgroundColor
Definition: Style.h:35
Float_t w
Definition: plot.C:20
auto zero()
Definition: PMNS.cxx:47
void ProjectionX(T *from, U *to)
Helper for Unweighted.
void next()
Definition: show_event.C:84
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:70
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
void drawBFMirror(double bfSin, double bfDm, Color_t color, Style_t marker, double size=1.5)
Draw best fit at both octants. Truth for neutrinos, not sure about antineutrinos. ...
Definition: Plots.cxx:1297
TPad * pad1
Definition: analysis.C:13
void GetSystBands(osc::IOscCalculator *calc, IPrediction *pred, std::vector< const ISyst * > systs, std::vector< TH1 * > &hUps, std::vector< TH1 * > &hDns, double pot)
Definition: Plots.cxx:245