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