Plots.cxx

Go to the documentation of this file.

#include "CAFAna/Analysis/Plots.h"

#include "CAFAna/Analysis/Style.h"

#include "CAFAna/Core/Ratio.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SystShifts.h"
#include "Utilities/func/MathUtil.h"

#include "TCanvas.h"
#include "TFeldmanCousins.h"
#include "TGraph.h"
#include "TGraphAsymmErrors.h"
#include "TH1.h"
#include "TH2.h"
#include "THStack.h"
#include "THLimitsFinder.h"
#include "TLegend.h"
#include "TLine.h"
#include "TList.h"
#include "TMarker.h"
#include "TPaveText.h"
#include "TText.h"

#include <algorithm>
#include <iostream>
#include <cmath>
#include <ctime>
#include <functional>
#include <memory>

namespace ana
{
  //----------------------------------------------------------------------
  TH1* DataMCComparison(const Spectrum& data, const Spectrum& mc, EBinType bintype)
  {
    TH1* ret = 0;

    TH1* hMC = mc.ToTH1(data.POT(), kTotalMCColor, kSolid, kPOT, bintype);

    TH1* hData = data.ToTH1(data.POT(), kBlack, kSolid, kPOT, bintype);
    hData->Sumw2();
    hData->SetMarkerStyle(kFullCircle);

    // Need to draw the biggest thing first to get correct axis limits
    if(hMC->GetMaximum() > hData->GetMaximum()+sqrt(hData->GetMaximum())){
      ret = hMC;
      hMC->Draw("hist");
      hData->Draw("ep same");
    }
    else{
      ret = hData;
      hData->Draw("ep");
      hMC->Draw("hist same");
      hData->Draw("ep same"); // Data always goes on top
    }

    gPad->Update();

    return ret;
  }

  //----------------------------------------------------------------------
  TH1* DataMCComparison(const Spectrum& data,
                   const IPrediction* mc,
                   osc::IOscCalc* calc,
                   const SystShifts & shifts,
                   EBinType bintype)
  {
    return DataMCComparison(data, shifts.IsNominal() ? mc->Predict(calc) : mc->PredictSyst(calc, shifts), bintype);
  }

  //----------------------------------------------------------------------
  TH1* DataMCComparisonAreaNormalized(const Spectrum& data, const Spectrum& mc)
  {
    TH1* ret = 0;

    TH1* hMC = mc.ToTH1(mc.POT());
    hMC->SetLineColor(kTotalMCColor);

    TH1* hData = data.ToTH1(data.POT());
    hData->Sumw2();
    hData->SetMarkerStyle(kFullCircle);

    hMC->Scale(hData->Integral()/hMC->Integral());

    // Need to draw the biggest thing first to get correct axis limits
    if(hMC->GetMaximum() > hData->GetMaximum()+sqrt(hData->GetMaximum())){
      ret = hMC;
      hMC->Draw("hist");
      hData->Draw("ep same");
    }
    else{
      ret = hData;
      hData->Draw("ep");
      hMC->Draw("hist same");
      hData->Draw("ep same"); // Data always goes on top
    }

    gPad->Update();

    return ret;
  }

  //----------------------------------------------------------------------
  TH1* DataMCComparisonAreaNormalized(const Spectrum& data,
                   const IPrediction* mc,
                   osc::IOscCalc* calc)
  {
    return DataMCComparisonAreaNormalized(data, mc->Predict(calc));
  }

  //----------------------------------------------------------------------
  TH1* DataMCComparisonComponents(const Spectrum& data,
                                  const IPrediction* mc,
                                  osc::IOscCalc* calc)
  {
    TH1* ret = 0;

    TH1* hMC = mc->Predict(calc).ToTH1(data.POT());
    hMC->SetLineColor(kTotalMCColor);

    TH1* hNC = mc->PredictComponent(calc,
                                    Flavors::kAll,
                                    Current::kNC,
                                    Sign::kBoth).ToTH1(data.POT());
    hNC->SetLineColor(kNCBackgroundColor);

    TH1* hCC = mc->PredictComponent(calc,
                                    Flavors::kAllNuMu,
                                    Current::kCC,
                                    Sign::kBoth).ToTH1(data.POT());
    hCC->SetLineColor(kNumuBackgroundColor);

    TH1* hBeam = mc->PredictComponent(calc,
                                      Flavors::kNuEToNuE,
                                      Current::kCC,
                                      Sign::kBoth).ToTH1(data.POT());
    hBeam->SetLineColor(kBeamNueBackgroundColor);

    TH1* hData = data.ToTH1(data.POT());
    hData->SetMarkerStyle(kFullCircle);

    // Need to draw the biggest thing first to get correct axis limits
    if(hMC->GetMaximum() > hData->GetMaximum()+sqrt(hData->GetMaximum())){
      ret = hMC;
      hMC->Draw("hist");
      hData->Draw("ep same");
    }
    else{
      ret = hData;
      hData->Draw("ep");
      hMC->Draw("hist same");
      hData->Draw("ep same"); // Data always goes on top
    }

    hNC->Draw("hist same");
    hCC->Draw("hist same");
    hBeam->Draw("hist same");

    gPad->Update();

    return ret;
  }

  //----------------------------------------------------------------------
  void DataMCRatio(const Spectrum& data,
                   const IPrediction* mc,
                   osc::IOscCalc* calc,
                   double miny, double maxy)
  {
    DataMCRatio(data, mc->Predict(calc), miny, maxy);
  }

  //----------------------------------------------------------------------
  void DataMCRatio(const Spectrum& data,
                   const Spectrum& mc,
                   double miny, double maxy)
  {
    Ratio ratio(data, mc);
    TH1* h = ratio.ToTH1();
    h->GetYaxis()->SetTitle("Data / MC");
    h->GetYaxis()->SetRangeUser(miny, maxy);
    h->Draw();

    TLine* one = new TLine(h->GetXaxis()->GetXmin(), 1,
                           h->GetXaxis()->GetXmax(), 1);
    one->SetLineWidth(2);
    one->SetLineStyle(7);
    one->Draw();

    gPad->Update();
  }

  //----------------------------------------------------------------------
  void DataMCAreaNormalizedRatio(const Spectrum& data,
                                 const IPrediction* mc,
                                 osc::IOscCalc* calc,
                                 double miny, double maxy)
  {
    DataMCAreaNormalizedRatio(data, mc->Predict(calc), miny, maxy);
  }

  //----------------------------------------------------------------------
  void DataMCAreaNormalizedRatio(const Spectrum& data,
                                 const Spectrum& mc,
                                 double miny, double maxy)
  {
    Spectrum mcScaled = mc;
    mcScaled.OverridePOT(mcScaled.POT() * mc.Integral(1) / data.Integral(1));

    DataMCRatio(data, mcScaled, miny, maxy);
  }

  //----------------------------------------------------------------------
  void RatioPlot(const Spectrum& data,
                 const Spectrum& expected,
                 const Spectrum& fit,
                 double miny, double maxy)
  {
    Ratio fitRatio(fit, expected);
    Ratio dataRatio(data, expected);

    TH1* h = fitRatio.ToTH1();
    h->GetYaxis()->SetTitle("Ratio to expectation");
    h->GetYaxis()->SetRangeUser(miny, maxy);
    h->SetLineColor(kTotalMCColor);
    h->Draw("][");

    h = dataRatio.ToTH1();
    h->SetMarkerStyle(kFullCircle);
    h->Draw("ep same");

    TLine* one = new TLine(h->GetXaxis()->GetXmin(), 1,
                           h->GetXaxis()->GetXmax(), 1);
    one->SetLineWidth(2);
    one->SetLineStyle(7);
    one->Draw();

    gPad->Update();
  }

  //----------------------------------------------------------------------
  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)
  {
    for (const ISyst* syst: systs) {
      SystShifts shifts;
    
      shifts.SetShift(syst, +1);
      TH1* hUp = pred->PredictComponentSyst(calc, shifts, flav, curr, sign).ToTH1(pot);
    
      shifts.SetShift(syst, -1);
      TH1* hDn = pred->PredictComponentSyst(calc, shifts, flav, curr, sign).ToTH1(pot);

      hUps.push_back(hUp);
      hDns.push_back(hDn);
    }

    return;
  }

  //----------------------------------------------------------------------
  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)
  {
    for (const ISyst* syst: systs) {
      SystShifts shifts;
    
      shifts.SetShift(syst, +1);
      TH1* hUp = pred->PredictComponentSyst(calc, shifts, flav, curr, sign).ToTH1(pot);
      hUp->Add( pred->PredictComponentSyst(calc, SystShifts::Nominal(), flav, curr, sign).ToTH1(pot), -1 ); // remove the nominal pred
      hUp->Add( pred->PredictComponentSyst(calc, bfshifts, flav, curr, sign).ToTH1(pot), +1 ); // add on the BF shifted pred
     
      shifts.SetShift(syst, -1);
      TH1* hDn = pred->PredictComponentSyst(calc, shifts, flav, curr, sign).ToTH1(pot);
      hDn->Add( pred->PredictComponentSyst(calc, SystShifts::Nominal(), flav, curr, sign).ToTH1(pot), -1 ); // remove the nominal pred
      hDn->Add( pred->PredictComponentSyst(calc, bfshifts, flav, curr, sign).ToTH1(pot), +1 ); // add on the BF shifted pred

      hUps.push_back(hUp);
      hDns.push_back(hDn);
    }

    return;
  }

  
  //----------------------------------------------------------------------
  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)
  {

    Spectrum nom = pred->Predict(calc);

    std::vector<Spectrum> ups, dns;

    for(const ISyst* syst: systs){
      SystShifts shifts;
      shifts.SetShift(syst, +1);
      ups.push_back(pred->PredictSyst(calc, shifts));
      shifts.SetShift(syst, -1);
      dns.push_back(pred->PredictSyst(calc, shifts));
    }

    return PlotWithSystErrorBand(nom, ups, dns, pot, col, errCol, headroom, newaxis, bintype);


  }

  //----------------------------------------------------------------------

  TGraphAsymmErrors*  PlotWithSystErrorBand(const Spectrum& nominal,
                                            const std::vector<Spectrum>& upShifts,
                                            const std::vector<Spectrum>& downShifts,
                                            double pot, int col, int errCol,
                                            float headroom, bool newaxis,
                                            EBinType bintype)
  {

    TH1* nom =  nominal.ToTH1(pot, kBlack, kSolid, kPOT, bintype);

    std::vector<TH1*> ups, dns;

    for(const auto& upShift:upShifts)     ups.push_back(upShift.ToTH1(pot, kBlack, kSolid, kPOT, bintype));
    for(const auto& downShift:downShifts) dns.push_back(downShift.ToTH1(pot, kBlack, kSolid, kPOT, bintype));

    TGraphAsymmErrors* g = PlotWithSystErrorBand(nom, ups, dns, col, errCol, headroom, newaxis);

    for(TH1* up: ups) delete up;
    for(TH1* dn: dns) delete dn;

    return g;
  }

  //----------------------------------------------------------------------

  TGraphAsymmErrors* PlotWithSystErrorBand(TH1*& nom,
                                           std::vector<TH1*>& ups,
                                           std::vector<TH1*>& dns,
                                           int col, int errCol,
                                           float headroom, bool newaxis)
  {
    if(col == -1){
      col = kTotalMCColor;
      errCol = kTotalMCErrorBandColor;
    }
    else if(errCol == -1) errCol = col-7; // hopefully a lighter version

    nom->SetLineColor(col);
    nom->GetXaxis()->CenterTitle();
    nom->GetYaxis()->CenterTitle();
    if(newaxis) nom->Draw("hist ]["); // Set the axes up

    double yMax = nom->GetBinContent(nom->GetMaximumBin());

    TGraphAsymmErrors* g = new TGraphAsymmErrors;

    for(int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
      const double y = nom->GetBinContent(binIdx);
      g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y);

      const double w = nom->GetXaxis()->GetBinWidth(binIdx);

      double errUp = 0;
      double errDn = 0;

      for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
        double hi = ups[systIdx]->GetBinContent(binIdx)-y;
        double lo = dns[systIdx]->GetBinContent(binIdx)-y;

        // If both shifts are in the same direction use the larger
        double min = std::min(hi,lo);
        double max = std::max(hi,lo);
        if(max < 0) max=0;
        if(min > 0) min=0;

        errUp += max*max;
        errDn += min*min;
      } // end for systIdx

      g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
    } // end for i

    g->SetFillColor(errCol);
    g->Draw("e2 same");
    g->GetYaxis()->SetRangeUser(0.00001, headroom*yMax);
    nom->GetYaxis()->SetRangeUser(0.0001, headroom*yMax);

    nom->Draw("hist same");

    gPad->RedrawAxis();
    
    return g;
  }

  //----------------------------------------------------------------------
  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)
  {

    Spectrum nom = pred->Predict(calc);

    std::vector<Spectrum> ups, dns;

    for(const ISyst* syst: systs){
      SystShifts shifts;
      shifts.SetShift(syst, +1);
      ups.push_back(pred->PredictSyst(calc, shifts));
      shifts.SetShift(syst, -1);
      dns.push_back(pred->PredictSyst(calc, shifts));
    }

    PlotWithAreaSystErrorBand(nom, ups, dns, pot, col, errCol, headroom, newaxis, bintype);


  }

  //----------------------------------------------------------------------
  void PlotWithAreaSystErrorBand(const Spectrum& nominal,
                             std::vector<Spectrum> upShifts,
                             std::vector<Spectrum> downShifts,
                             double pot, int col, int errCol,
                             float headroom, bool newaxis,
                                            EBinType bintype)
  {
    double nomIntegral = nominal.Integral(pot);
    if (nomIntegral > 0)
    {
      // ugh so much boilerplate
      for (auto & collection : std::vector<std::reference_wrapper<std::vector<Spectrum>>>{ upShifts, downShifts })
      {
        for (Spectrum & shiftSpec : collection.get())
        {
          // this is a dirty trick, but these Spectra are going to be thrown out once this method exits anyway
          shiftSpec.OverridePOT(shiftSpec.POT() * shiftSpec.Integral(pot) / nomIntegral);
        }
      }
    }

    PlotWithSystErrorBand(nominal, upShifts, downShifts, pot, col, errCol, headroom, newaxis, bintype);
  }

  //----------------------------------------------------------------------
  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)
  {
    TH1* nom =  nominal.ToTH1(pot, kBlack, kSolid, kPOT, bintype);
    TH1* alt =  alternative.ToTH1(pot, kBlack, kSolid, kPOT, bintype);
    std::vector<TH1*> ups, dns;
    std::vector<TH1*> upsA, dnsA;
    for(const auto& upShift:upShifts)     ups.push_back(upShift.ToTH1(pot, kBlack, kSolid, kPOT, bintype));
    for(const auto& downShift:downShifts) dns.push_back(downShift.ToTH1(pot, kBlack, kSolid, kPOT, bintype));
    for(const auto& upShiftAlt:upShiftsAlt)     upsA.push_back(upShiftAlt.ToTH1(pot, kBlack, kSolid, kPOT, bintype));
    for(const auto& downShiftAlt:downShiftsAlt) dnsA.push_back(downShiftAlt.ToTH1(pot, kBlack, kSolid, kPOT, bintype));
    nom->SetLineColor(col1);
    nom->GetXaxis()->CenterTitle();
    nom->GetYaxis()->CenterTitle();
    if(newaxis) nom->Draw("hist"); // Set the axes up
    alt->SetLineColor(col2);
    double yMax = nom->GetBinContent(nom->GetMaximumBin());
    TGraphAsymmErrors* g = new TGraphAsymmErrors;
    TGraphAsymmErrors* g2 = new TGraphAsymmErrors;
    for(int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
      const double y = nom->GetBinContent(binIdx);
      const double y2 = alt->GetBinContent(binIdx);
      g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y);
      g2->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y2);
      const double w = nom->GetXaxis()->GetBinWidth(binIdx);
      double errUp = 0, errUp2 = 0;
      double errDn = 0, errDn2 = 0;
      for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
        double hi = ups[systIdx]->GetBinContent(binIdx)-y;
        double hi2 = upsA[systIdx]->GetBinContent(binIdx)-y2;
        double lo = y-dns[systIdx]->GetBinContent(binIdx);
        double lo2 = y2-dnsA[systIdx]->GetBinContent(binIdx);
        if(lo <= 0 && hi <= 0) std::swap(lo, hi);
        if(lo2 <= 0 && hi2 <= 0) std::swap(lo2, hi2);
        errUp += hi*hi;
        errDn += lo*lo;
        errUp2 += hi2*hi2;
        errDn2 += lo2*lo2;
        // TODO: what happens if they're both high or both low?
      } // end for systIdx
      g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
      g2->SetPointError(binIdx, w/2, w/2, sqrt(errDn2), sqrt(errUp2));
    } // end for i
    g2->SetFillColorAlpha(col2, 0.3);
    g2->Draw("2");
    // g2->GetYaxis()->SetRangeUser(0, headroom*yMax);
    g2->GetYaxis()->SetRangeUser(0.0, 19.0);
    g->SetFillColorAlpha(col1, 0.3);
    g->Draw("2");
    // g->GetYaxis()->SetRangeUser(0, headroom*yMax);
    // nom->GetYaxis()->SetRangeUser(0, headroom*yMax);
    g->GetYaxis()->SetRangeUser(0.0, yMax*headroom);
    nom->GetYaxis()->SetRangeUser(0.0, yMax*headroom);
    alt->Draw("hist same");
    nom->Draw("hist same");
    for(TH1* up: ups) delete up;
    for(TH1* dn: dns) delete dn;
    for(TH1* up2: upsA) delete up2;
    for(TH1* dn2: dnsA) delete dn2;
  }





  //----------------------------------------------------------------------
  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)
  {

    Spectrum nom = pred->Predict(calc);

    std::vector<Spectrum> ups, dns;

    for(const ISyst* syst: systs){
      SystShifts shifts;
      shifts.SetShift(syst, +1);
      ups.push_back(pred->PredictSyst(calc, shifts));
      shifts.SetShift(syst, -1);
      dns.push_back(pred->PredictSyst(calc, shifts));
    }

    return PlotWithSystErrorBand_Quant(quant, nom, ups, dns, pot, col, errCol, maxy, newaxis);


  }

  //----------------------------------------------------------------------
  TGraphAsymmErrors*  PlotWithSystErrorBand_Quant(const int quant,
                                                  const Spectrum& nominal,
                                                  const std::vector<Spectrum>& upShifts,
                                                  const std::vector<Spectrum>& downShifts,
                                                  double pot, int col, int errCol,
                                                  float maxy, bool newaxis)
  {

    TH1* nom =  nominal.ToTH1(pot);
    nom->GetYaxis()->SetTitle("Events/0.1 GeV");
    std::vector<TH1*> ups, dns;

    for(const auto& upShift:upShifts)     ups.push_back(upShift.ToTH1(pot));
    for(const auto& downShift:downShifts) dns.push_back(downShift.ToTH1(pot));
    
    TGraphAsymmErrors* g = PlotWithSystErrorBand_Quant(quant, nom, ups, dns, col, errCol, maxy, newaxis);
    
    for(TH1* up: ups) delete up;
    for(TH1* dn: dns) delete dn;

    return g;
  }

  //----------------------------------------------------------------------
  TGraphAsymmErrors* PlotWithSystErrorBand_Quant(const int quant,
                                                 TH1*& nom,
                                                 std::vector<TH1*>& ups,
                                                 std::vector<TH1*>& dns,
                                                 int col, int errCol,
                                                 float maxy, bool newaxis)
  {
    if(col == -1){
      col = kTotalMCColor;
      errCol = kTotalMCErrorBandColor;
    }
    else if(errCol == -1) errCol = col-7; // hopefully a lighter version

    nom->GetYaxis()->SetTitle("Events/0.1 GeV");
    nom->SetLineColor(col);
    nom->GetXaxis()->CenterTitle();
    nom->GetYaxis()->CenterTitle();
    if(newaxis) nom->Draw("hist ]["); // Set the axes up

    TGraphAsymmErrors* g = new TGraphAsymmErrors;

    for(int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
      const double y = nom->GetBinContent(binIdx);
      g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y);

      const double w = nom->GetXaxis()->GetBinWidth(binIdx);

      double errUp = 0;
      double errDn = 0;

      for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
        double hi = ups[systIdx]->GetBinContent(binIdx)-y;
        double lo = dns[systIdx]->GetBinContent(binIdx)-y;

        // If both shifts are in the same direction use the larger
        double min = std::min(hi,lo);
        double max = std::max(hi,lo);
        if(max < 0) max=0;
        if(min > 0) min=0;

        errUp += max*max;
        errDn += min*min;
      } // end for systIdx

      g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
    } // end for i


    // setting 2x2 = 4 pads canvas
    if(quant==1 || quant==2){
      nom->GetXaxis()->SetLabelSize(0.);
    }
    if(quant==2 || quant==4){
      nom->GetYaxis()->SetLabelSize(0.);
    }
    float minx = 0.99 * nom->GetXaxis()->GetXmin();
    float maxx = 0.99 * nom->GetXaxis()->GetXmax();
    nom->GetXaxis()->SetLabelOffset(0.01);
    nom->GetYaxis()->SetLabelOffset(0.01);
    nom->GetYaxis()->SetTitleOffset(0.50);
    nom->GetXaxis()->SetTitleSize(0.0);
    nom->GetYaxis()->SetTitleSize(0.0);
    nom->GetXaxis()->SetTitleSize(0.0);
    nom->GetYaxis()->SetTitleSize(0.0);
    nom->GetYaxis()->SetRangeUser(0.00001, maxy);
    nom->GetXaxis()->SetRangeUser(minx, maxx);
    g->GetYaxis()  ->SetRangeUser(0.00001, maxy);
    g->GetXaxis()  ->SetRangeUser(minx, maxx);
    g->SetFillColor(errCol);
    g->Draw("e2 same");
    nom->Draw("hist ][ same");

    return g;
  }






  //----------------------------------------------------------------------
  THStack* ToTHStack(const std::vector<std::pair<const Spectrum&, Color_t>>& s,
                     double pot)
  {
    THStack* ret = new THStack;
    for(auto it: s){
      TH1* h = it.first.ToTH1(pot);
      h->SetFillColor(it.second);
      ret->Add(h, "hist");
    }
    return ret;
  }

  //----------------------------------------------------------------------
  /// Helper for \ref AutoPlaceLegend
  double PointDistanceToBox(double x, double y,
                            double x0, double y0, double x1, double y1)
  {
    // Inside
    if(x > x0 && x < x1 && y > y0 && y < y1) return 0;

    // Corners
    double d = util::sqr(x-x0)+util::sqr(y-y0);
    d = std::min(d, util::sqr(x-x1)+util::sqr(y-y0));
    d = std::min(d, util::sqr(x-x1)+util::sqr(y-y1));
    d = std::min(d, util::sqr(x-x0)+util::sqr(y-y1));

    // Top and bottom edges
    if(x > x0 && x < x1){
      d = std::min(d, util::sqr(y-y0));
      d = std::min(d, util::sqr(y-y1));
    }
    // Left and right
    if(y > y0 && y < y1){
      d = std::min(d, util::sqr(x-x0));
      d = std::min(d, util::sqr(x-x1));
    }

    return d;
  }

  //----------------------------------------------------------------------
  TLegend* AutoPlaceLegend(double dx, double dy, double yPin)
  {
    gPad->Update();

    // Convert requested width and height into physics coordinates
    dx *= (gPad->GetX2()-gPad->GetX1());
    dy *= (gPad->GetY2()-gPad->GetY1());

    // Range of axes in physics units
    const double x0 = gPad->GetUxmin();
    const double x1 = gPad->GetUxmax();
    const double y0 = gPad->GetUymin();
    const double y1 = gPad->GetUymax();

    const double X = x1-x0;
    const double Y = y1-y0;

    double bestd = 0;
    double bestx = 0;
    double besty = 0;
    // If we want to pin Y pos, set it to that now.
    if(yPin >= 0) besty = yPin * (gPad->GetY2()-gPad->GetY1());

    for(bool fallback: {false, true}){
      for(double x = x0+dx/2; x < x1-dx/2; x += X/50){
        for(double y = y0+dy/2; y < y1-dy/2; y += Y/50){
          double d = 999999;

          // Repel from edges
          d = std::min(d, util::sqr((x-dx/2-x0)/X));
          d = std::min(d, util::sqr((x+dx/2-x1)/X));
          d = std::min(d, util::sqr((y-dy/2-y0)/Y));
          d = std::min(d, util::sqr((y+dy/2-y1)/Y));
          if(d < bestd) continue;

          TIter next(gPad->GetListOfPrimitives());
          while(TObject* obj = next()){

            if(!obj->InheritsFrom(TH1::Class())) continue;
            if( obj->InheritsFrom(TH2::Class())) continue;

            TH1* h = (TH1*)obj;

            for(int n = 1; n <= h->GetNbinsX(); ++n){
              const double px = h->GetBinCenter(n);
              const double py = h->GetBinContent(n);

              if(fallback){
                d = std::min(d, util::sqr((px-x)/X)+util::sqr((py-y)/Y));
              }
              else{
                d = std::min(d, PointDistanceToBox(px/X, py/Y,
                                                   (x-dx/2)/X, (y-dy/2)/Y,
                                                   (x+dx/2)/X, (y+dy/2)/Y));
              }
              if(d < bestd) break;
            }
            if(d < bestd) break;
          } // end while

          if(d > bestd){
            bestd = d;
            bestx = x;
            // Update Y if we're not pinning it.
            if (yPin < 0) besty = y;
          }
        } // end for y
      } // end for x

      if(bestd != 0) break; // If we always collide, have to do fallback
    } // end for fallback

    // Convert to pad coordinates
    const double nx = (bestx-gPad->GetX1())/(gPad->GetX2()-gPad->GetX1());
    const double ny = (besty-gPad->GetY1())/(gPad->GetY2()-gPad->GetY1());

    const double ndx = dx/(gPad->GetX2()-gPad->GetX1());
    const double ndy = dy/(gPad->GetY2()-gPad->GetY1());

    return new TLegend(nx-ndx/2, ny-ndy/2, nx+ndx/2, ny+ndy/2);
  }

  //----------------------------------------------------------------------
  void ErrorBarChart(const std::map<std::string, std::pair<double, double>>& systErrors,
                     const std::pair<double, double> & statErr,
                     const std::string & label )
  {
    bool doStat = statErr.first != 0 && statErr.second != 0;
    int nbins = systErrors.size() + 1;  // extra for the "total systematic" bar

    // Systematics are on top, total systematic and (optional) statistical
    // error bar are at the bottom.
    int firstSystBin = 1;
    if(doStat){
      ++nbins;
      ++firstSystBin;
    }

    // Sum all the systematics in qaudrature
    std::pair<double, double> systTot = std::make_pair(0., 0.);
    for(auto it_systPairs: systErrors)
    {
      systTot.first += util::sqr(it_systPairs.second.first);
      systTot.second += util::sqr(it_systPairs.second.second);
    }
    systTot.first = sqrt(systTot.first);
    systTot.second = sqrt(systTot.second);

    // x-axis range is a multiple of 5% with at least 2% padding
    double xrange = 0;
    while(xrange < std::max(std::max(systTot.first, systTot.second), std::max(statErr.first, statErr.second))+2) xrange += 5;

    TH2* axes = new TH2F("", (";" + label).c_str(),
                         10, -xrange, +xrange, nbins, 0, nbins);
    axes->GetXaxis()->CenterTitle();

    // Put the systematics in a vector so we can sort them by size
    std::vector<std::pair<std::pair<double, double>, std::string>> systList;
    for(auto it: systErrors) systList.push_back(std::make_pair(it.second, it.first));
    std::sort(systList.begin(),
              systList.end(),
              [](const auto & errNamePair1, const auto & errNamePair2)
              {
                return ( std::max(errNamePair1.first.first, errNamePair1.first.second)
                         < std::max(errNamePair2.first.first, errNamePair2.first.second) );
              }
    );

    // Label the y axis
    TAxis* yax = axes->GetYaxis();
    for(unsigned int i = 0; i < systList.size(); ++i)
      yax->SetBinLabel(i+firstSystBin+1, systList[i].second.c_str());

    yax->SetBinLabel(firstSystBin, "Total syst. error");
    if(doStat) yax->SetBinLabel(1, "Statistical error");

    // Make the labels approximately match the x axis title
    yax->SetLabelSize(.06);

    axes->Draw();

    // Blue boxes for each systematic
    for(unsigned int i = 0; i < systList.size(); ++i){
      TBox* box = new TBox(-systList[i].first.first,  i+firstSystBin+.2,
                           +systList[i].first.second, i+firstSystBin+.8);
      std::cout<<"ERROR = (" << systList[i].first.first << ", " << systList[i].first.second << ")" << std::endl;
      box->SetFillColor(kAzure-8);
      box->Draw();
    }

    // Total systematic error
    TBox* syst = new TBox(-systTot.first, firstSystBin-.8,
                          +systTot.second, firstSystBin-.2);
    std::cout<<"TOTAL SYSTEMATIC ERROR = (" << systTot.first << ", " << systTot.second << ")" << std::endl;
    syst->SetFillColor(kAzure-8);
    syst->Draw();

    // Red box for statistical error
    if(doStat){
      TBox* stat = new TBox(-statErr.first, .2, +statErr.second, .8);
      stat->SetFillColor(kRed-7);
      stat->Draw();
    }

    // Separate individual systematics from the total and statistical errors
    TLine* div = new TLine(-xrange, firstSystBin, +xrange, firstSystBin);
    div->SetLineWidth(2);
    div->Draw();

    // Vertical line marking the symmetry point
    TLine* zero = new TLine(0, 0, 0, nbins);
    zero->SetLineStyle(7);
    zero->SetLineWidth(2);
    zero->Draw();

    // Leave enough space for the systematic labels
    gPad->SetLeftMargin(.25);

  }

  //----------------------------------------------------------------------
  void CountingExperimentErrorBarChart(const std::map<std::string, double>& systbars, double statErr, bool bkgdOrSig, bool shortchart)
  {
    std::map<std::string, double> systs = systbars; // breaks shortchart option (shortchart) ? SumNueSecondAnaSysts(systbars) : systbars;

    std::map<std::string, std::pair<double, double>> systVals;
    for (const auto & systPair : systbars)
      systVals[systPair.first] = std::make_pair(systPair.second, systPair.second);
    auto statVals = std::make_pair(statErr, statErr);

    std::string label = std::string(bkgdOrSig ? "Signal" : "Background") + " uncertainty (%)";
    ErrorBarChart(systVals, statVals, label);
  }

  //----------------------------------------------------------------------
  TGraphAsymmErrors* GraphWithPoissonErrors(const TH1* h, bool noErrorsXaxis, bool drawEmptyBins)
  {
    TGraphAsymmErrors* gr = new TGraphAsymmErrors(h);

    TFeldmanCousins fc(0.6827);//1 sigma

    for(int i = 0; i < h->GetNbinsX(); ++i){
      double x, y;
      gr->GetPoint(i, x, y);

      if ( drawEmptyBins || y!=0 ){
        if ( y < 50 )   gr->SetPointEYlow(i, y-fc.CalculateLowerLimit(y,0));
        else            gr->SetPointEYlow(i,sqrt(y));
        if ( y < 30 )   gr->SetPointEYhigh(i,fc.CalculateUpperLimit(y,0)-y);
        else            gr->SetPointEYhigh(i,sqrt(y));
      }

      if(noErrorsXaxis){
        gr->SetPointEXlow(i, 0);
        gr->SetPointEXhigh(i, 0);
      } // Do not use bin width as X-error for points
    }
    gr->SetMarkerStyle(kFullCircle);
    gr->SetMarkerColor(1);
    gr->SetMarkerSize(1);
    gr->SetLineWidth(2);


    return gr;
  }

  //----------------------------------------------------------------------
  TGraphAsymmErrors* GraphWithPoissonErrors2(const TH1* h, const TH1* h2, bool noErrorsXaxis, bool drawEmptyBins)
  {
    TGraphAsymmErrors* gr = new TGraphAsymmErrors(h);
    TGraphAsymmErrors* gr2 = new TGraphAsymmErrors();

    TFeldmanCousins fc(0.6827);//1 sigma

    int j = 0;
    for(int i = 0; i < h->GetNbinsX(); ++i){
      double x, y;
      gr->GetPoint(i, x, y);

      y = round(y);

      if ( (drawEmptyBins && h2->GetBinContent(i+1) > 0.05) || y!=0 ){
        gr2->SetPoint(j, x, y);
        gr2->SetPointEXlow(j, gr->GetErrorXlow(i));
        gr2->SetPointEXhigh(j, gr->GetErrorXhigh(i));
        if ( y < 50 )   gr2->SetPointEYlow(j, y-fc.CalculateLowerLimit(y,0));
        else            gr2->SetPointEYlow(j,sqrt(y));
        if ( y < 30 )   gr2->SetPointEYhigh(j,fc.CalculateUpperLimit(y,0)-y);
        else            gr2->SetPointEYhigh(j,sqrt(y));
        j++;
      }

      if(noErrorsXaxis){
        gr2->SetPointEXlow(i, 0);
        gr2->SetPointEXhigh(i, 0);
      } // Do not use bin width as X-error for points
    }

    return gr2;
  }


  //----------------------------------------------------------------------
  TGraph* graphAsymmErrorScaled(TH1* histScaled, TH1* hist, double overallScale){

    for(int bin = 1; bin <= hist->GetNbinsX(); bin++) {
      if(hist->GetBinContent(bin)==0){
        hist->SetBinContent(bin,0.001);
        histScaled->SetBinContent(bin,0.001);
      }
    }

    TFeldmanCousins fc(0.6827); // to get correct poisson error bars
    TGraphAsymmErrors *datapoisson = new TGraphAsymmErrors(histScaled);

    for(int bin = 1; bin <= hist->GetNbinsX(); bin++) {
      //default numu: scale by bin 0.1/width
      double binWidth = hist->GetBinWidth(bin);
      double scaleFactor = overallScale / binWidth;
      double y = hist->GetBinContent(bin);

      // since this is data it should be a whole number
      // fc Upper/Lower limit is finicky eg if y = 3.01
      // the returned value is the same as when y = 4
      y = round(y);

      double errup = 0, errdn = 0;
      if(y < 30) errup = (fc.CalculateUpperLimit(y,0) - y);
      else       errup = sqrt(y);
      if(y < 50) errdn = (y - fc.CalculateLowerLimit(y,0));
      else       errdn = sqrt(y);

      datapoisson->SetPointEYhigh(bin-1, scaleFactor * errup);
      datapoisson->SetPointEYlow (bin-1, scaleFactor * errdn);
    }
    datapoisson->SetMarkerStyle(kFullCircle);
    datapoisson->SetMarkerColor(1);
    datapoisson->SetMarkerSize(1);
    datapoisson->SetLineWidth(2);

    return datapoisson;

  }

  //----------------------------------------------------------------------
  TGraph* graphAsymmErrorWithBkgScaled(TH1* data, TH1* bkgd, double overallScale){

    for(int bin = 1; bin <= data->GetNbinsX(); bin++) {
      if(data->GetBinContent(bin)==0){
        data->SetBinContent(bin,0.001);
      }
    }

    TFeldmanCousins fc(0.6827); // to get correct poisson error bars
    TGraphAsymmErrors *gr = new TGraphAsymmErrors(data);

    for(int bin = 1; bin <= data->GetNbinsX(); bin++) {
      //default numu: scale by bin 0.1/width
      double width = data->GetBinWidth(bin);
      double scale = width / overallScale; // this is to reverse scaling

      double x, y;
      gr->GetPoint(bin-1, x, y);
      double bkg = bkgd->GetBinContent(bin);
      y   *= scale;
      bkg *= scale;

      // since this is data it should be a whole number
      // fc Upper/Lower limit is finicky eg if y = 3.01
      // the returned value is the same as when y = 4
      y = round(y);

      double errup = 0, errdn = 0;
      if(y < 30) errup = (fc.CalculateUpperLimit(y,bkg) - y);
      else       errup = sqrt(y);
      if(y < 50) errdn = (y - fc.CalculateLowerLimit(y,bkg));
      else       errdn = sqrt(y);

      gr->SetPointEYhigh(bin-1, 1/scale * errup);
      gr->SetPointEYlow (bin-1, 1/scale * errdn);
    }

    gr->SetMarkerStyle(kFullCircle);
    gr->SetMarkerColor(1);
    gr->SetMarkerSize(1);
    gr->SetLineWidth(2);

    return gr;

  }

  //----------------------------------------------------------------------
  TGraph* RatioAsymmError(TH1* data, TH1* pred)
  {
    TGraphAsymmErrors* gRat = GraphWithPoissonErrors(data);
    TFeldmanCousins fc(0.6827); // to get correct poisson error bars

    // Hist underflow bin is 0 so index from 1
    // Graphs start at 0 though so will offset points
    for (int bin = 1; bin <= data->GetNbinsX(); ++bin) {
      double x, y;
      gRat->GetPoint(bin-1, x, y);
      double den  = pred->GetBinContent(bin);

      // since this is data it should be a whole number
      // fc Upper/Lower limit is finicky eg if y = 3.01
      // the returned value is the same as when y = 4
      // Also negatives are weird
      y = round(y);
      if (y <= 0) y = 0.01;

      double errup = 0, errdn = 0;
      if (y < 30) errup = (fc.CalculateUpperLimit(y, 0) - y) / den;
      else        errup = sqrt(y) / den;
      if (y < 50) errdn = (y - fc.CalculateLowerLimit(y, 0)) / den;
      else        errdn = sqrt(y) / den;
      
      gRat->SetPoint(bin-1, x, y / den);
      gRat->SetPointEYhigh(bin-1, errup);
      gRat->SetPointEYlow(bin-1, errdn);
    }
    gRat->SetMarkerStyle(kFullCircle);
    gRat->SetMarkerColor(1);
    gRat->SetMarkerSize(1);
    gRat->SetLineWidth(2);
    
    return gRat;
  }

  //----------------------------------------------------------------------
  TGraph* RatioAsymmErrorScaled(TH1* data, TH1* pred, double overallScale)
  {
    TGraphAsymmErrors* gRat = GraphWithPoissonErrors(data);
    TFeldmanCousins fc(0.6827); // to get correct poisson error bars

    // Hist underflow bin is 0 so index from 1
    // Graphs start at 0 though so will offset points
    for (int bin = 1; bin <= data->GetNbinsX(); ++bin) {
      double width = data->GetBinWidth(bin);
      double scale = width / overallScale; // this is to reverse scaling

      double x, y;
      gRat->GetPoint(bin-1, x, y);
      double den  = pred->GetBinContent(bin);
      y   *= scale;
      den *= scale;

      // since this is data it should be a whole number
      // fc Upper/Lower limit is finicky eg if y = 3.01
      // the returned value is the same as when y = 4
      // Also negatives are weird
      y = round(y);
      if (y <= 0) y = 0.01;
      
      double errup = 0, errdn = 0;
      if (y < 30) errup = (fc.CalculateUpperLimit(y, 0) - y) / den;
      else        errup = sqrt(y) / den;
      if (y < 50) errdn = (y - fc.CalculateLowerLimit(y, 0)) / den;
      else        errdn = sqrt(y) / den;
      
      gRat->SetPoint(bin-1, x, y / den);
      gRat->SetPointEYhigh(bin-1, errup);
      gRat->SetPointEYlow (bin-1, errdn);
    }
    gRat->SetMarkerStyle(kFullCircle);
    gRat->SetMarkerColor(1);
    gRat->SetMarkerSize(1);
    gRat->SetLineWidth(2);
    
    return gRat;
  }

  //----------------------------------------------------------------------
  TGraph* RatioAsymmErrorWithBkg(TH1* data, TH1* pred, TH1* bkgd)
  {
    TGraphAsymmErrors* gRat = GraphWithPoissonErrors(data);
    TFeldmanCousins fc(0.6827); // to get correct poisson error bars

    // Hist underflow bin is 0 so index from 1
    // Graphs start at 0 though so will offset points
    for (int bin = 1; bin <= data->GetNbinsX(); ++bin) {
      double x, y;
      gRat->GetPoint(bin-1, x, y);
      double den = pred->GetBinContent(bin);
      double bkg = bkgd->GetBinContent(bin);

      // since this is data it should be a whole number
      // fc Upper/Lower limit is finicky eg if y = 3.01
      // the returned value is the same as when y = 4
      // Also negatives are weird
      y = round(y);
      if (y <= 0) y = 0.01;
      
      double errup = 0, errdn = 0;
      if (y < 30) errup = (fc.CalculateUpperLimit(y, bkg) - y) / den;
      else        errup = sqrt(y) / den;
      if (y < 50) errdn = (y - fc.CalculateLowerLimit(y, bkg)) / den;
      else        errdn = sqrt(y) / den;
      
      gRat->SetPoint(bin-1, x, y / den);
      gRat->SetPointEYhigh(bin-1, errup);
      gRat->SetPointEYlow (bin-1, errdn);
    }
    gRat->SetMarkerStyle(kFullCircle);
    gRat->SetMarkerColor(1);
    gRat->SetMarkerSize(1);
    gRat->SetLineWidth(2);
    
    return gRat;
  }

  
  //----------------------------------------------------------------------
  TGraph* RatioAsymmErrorWithBkgScaled(TH1* data, TH1* pred, TH1* bkgd, double overallScale)
  {
    TGraphAsymmErrors* gRat = GraphWithPoissonErrors(data); 
    TFeldmanCousins fc(0.6827); // to get correct poisson error bars

    // Hist underflow bin is 0 so index from 1
    // Graphs start at 0 though so will offset points
    for (int bin = 1; bin <= data->GetNbinsX(); ++bin) {
      double width = data->GetBinWidth(bin);
      double scale = width / overallScale; // this is to reverse scaling

      double x, y;
      gRat->GetPoint(bin-1, x, y);
      double den = pred->GetBinContent(bin);
      double bkg = bkgd->GetBinContent(bin);
      
      y   *= scale;
      den *= scale;
      bkg *= scale;

      // since this is data it should be a whole number
      // fc Upper/Lower limit is finicky eg if y = 3.01
      // the returned value is the same as when y = 4
      // Also negatives are weird
      y = round(y);
      if (y <= 0) y = 0.01;
            
      double errup = 0, errdn = 0;
      if (y < 30) errup = (fc.CalculateUpperLimit(y, bkg) - y) / den;
      else        errup = sqrt(y) / den;
      if (y < 50) errdn = (y - fc.CalculateLowerLimit(y, bkg)) / den;
      else        errdn = sqrt(y) / den;
      
      gRat->SetPoint(bin-1, x, y / den);
      gRat->SetPointEYhigh(bin-1, errup);
      gRat->SetPointEYlow (bin-1, errdn);
    }
    gRat->SetMarkerStyle(kFullCircle);
    gRat->SetMarkerColor(1);
    gRat->SetMarkerSize(1);
    gRat->SetLineWidth(2);
    
    return gRat;
  }

  //----------------------------------------------------------------------
  TGraph* ShadeBetweenHistograms(TH1* hmin, TH1* hmax)
  {
    int n = hmin->GetNbinsX();
    TGraph* gr = new TGraph(4*n);

    for(int i=1; i<=n; i++)
      {
        double xdown = hmax->GetBinLowEdge(i);
        double xup = hmax->GetBinLowEdge(i+1);
        double ydown = hmin->GetBinContent(i);
        double yup = hmax->GetBinContent(i);

        gr->SetPoint(2*(i-1), xdown, ydown);
        gr->SetPoint(2*(i-1)+1, xup, ydown);
        gr->SetPoint(4*n-2*i, xup, yup);
        gr->SetPoint(4*n-2*i+1, xdown, yup);
      }

      gr->SetFillColor(hmin->GetLineColor() - 9); // In principle this is lighter
      gr->SetLineColor(hmin->GetLineColor() - 9); // In case one does Draw("lf")

      return gr;
  }
  //----------------------------------------------------------------------
  TGraphAsymmErrors * ProfileQuantile(const TH2 * hist,
                                      const std::string & axisName,
                                      const std::string & graphName,
                                      const std::pair<double, double> & quantileDivisions)
  {
    if (hist->GetDimension() != 2)
            throw std::runtime_error(Form("Can't profile a histogram with other than 2 dimensions.  Yours is a %d-D histogram...", hist->GetDimension()));

    TAxis const* axis = nullptr;
    std::function<TH1D*(const char *, Int_t, Int_t, Option_t*)> projectionMethod;
    if (axisName == "x" || axisName == "X")
    {
      axis = hist->GetXaxis();
      projectionMethod = std::bind(&TH2::ProjectionY, hist, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4);
    }
    else if (axisName == "y" || axisName == "Y")
    {
      axis = hist->GetYaxis();
      projectionMethod = std::bind(&TH2::ProjectionX, hist, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4);
    }
    else
      throw std::runtime_error( Form("Can't profile onto unknown axis: '%s'", axisName.c_str()) );


    std::vector<double> points_x, points_y, errors_x_low, errors_x_high, errors_y_low, errors_y_high;

    double quantiles[2] = {0, 0};
    double _quantileDivisions[2] = {quantileDivisions.first, quantileDivisions.second};
    for (int binNum = 1; binNum <= axis->GetNbins(); binNum++)  // remember, 0 is underflow and Nbins+1 is overflow...
    {
      std::unique_ptr<TH1D> proj ( projectionMethod(Form("%s_tmp_%d", hist->GetName(), int(std::time(nullptr))), binNum, binNum, "") );  // randomish name, just to be safe

      double y = proj->GetMean();
      points_x.push_back(axis->GetBinCenter(binNum));
      points_y.push_back(y);

      // for now, make the errors for an empty projection 0
      unsigned int n_qs = 0;
      if (proj->Integral() == 0)
      {
        n_qs = 2;
        quantiles[0] = quantiles[1] = 0;
      }
      else
        n_qs = proj->GetQuantiles(2, quantiles, _quantileDivisions);
      if (n_qs != 2)
        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) );

      double binWidth = axis->GetBinWidth(binNum);
      errors_x_low.push_back(binWidth/2.);
      errors_x_high.push_back(binWidth/2.);
      errors_y_low.push_back(y-quantiles[0]);
      errors_y_high.push_back(quantiles[1]-y);
    }

    TGraphAsymmErrors * outGraph = new TGraphAsymmErrors(
            points_x.size(),
            &points_x[0],
            &points_y[0],
            &errors_x_low[0],
            &errors_x_high[0],
            &errors_y_low[0],
            &errors_y_high[0]
    );
    std::string name = (graphName.size()) ? graphName : Form("%s_quantile_errors", hist->GetName());
    outGraph->SetName(name.c_str());

    return outGraph;

  }


  // Draw single BF
  void drawBFSingle(double bfSin, double bfDm, Color_t color, Style_t marker, double size = 1.5)
  {

    TMarker* manMarker = new TMarker(bfSin, bfDm, marker);
    manMarker->SetMarkerSize(size);
    manMarker->SetMarkerColor(color);
    manMarker->Draw();

  }

  // Truth for FHC; not sure about RHC
  void drawBFMirror(double bfSin, double bfDm, Color_t color, Style_t marker, double size = 1.5)
  {

    double mirror = bfSin;
    TMarker* manMarker = new TMarker(bfSin, bfDm, marker);
    manMarker->SetMarkerSize(size);
    manMarker->SetMarkerColor(color);
    manMarker->Draw();

    if(bfSin > 0.514) mirror = 0.514 - (bfSin - 0.514);
    if(bfSin < 0.514) mirror = 0.514 + (0.514 - bfSin );
    TMarker* mirrorMarker = new TMarker(mirror, bfDm, marker);
    mirrorMarker->SetMarkerSize(size);
    mirrorMarker->SetMarkerColor(color);
    mirrorMarker->Draw();

  }



  //----------------------------------------------------------------------
  void MakeHistCanvasReady_Quant(const int quant,
                                 TH1* hist,
                                 double ymax){
    
    hist->GetXaxis()->SetLabelOffset(0.01);
    hist->GetYaxis()->SetLabelOffset(0.01);
    hist->GetYaxis()->SetTitleOffset(0.50);
    hist->GetXaxis()->SetTitleSize(0.0);
    hist->GetYaxis()->SetTitleSize(0.0);
    hist->GetYaxis()->SetRangeUser(0.0, ymax);

    if(quant==1 || quant==2) hist->GetXaxis()->SetLabelSize(0.);
    if(quant==2 || quant==4) hist->GetYaxis()->SetLabelSize(0.);
    // Q3 might need some axis labels removed to prevent interference.
    // Rather than assuming there will be trouble we can check the
    // optimization that ROOT does when plottng the axes
    // ChangeLabel labnum -1 alters the last label
    if(quant == 3) {
      // initialize the optimized min, max, step size,
      // and number of steps of the labels
      double optmin = 0., optmax = 0., optwi = 0.;
      int   optstp = 0;
                
      double xmin = hist->GetXaxis()->GetXmin();
      double xmax = hist->GetXaxis()->GetXmax();
      int   xdiv = hist->GetXaxis()->GetNdivisions() % 100;
      THLimitsFinder::Optimize(xmin, xmax, xdiv, optmin, optmax, optstp, optwi);

      if ( (xmax - optmax) < optwi/2 )
        hist->GetXaxis()->ChangeLabel( -1, -1, 0);

      optmin = 0.; optmax = 0.; optwi = 0.; optstp = 0;
      double ymin = 0.;
      int   ydiv = hist->GetYaxis()->GetNdivisions() % 100;
      THLimitsFinder::Optimize(ymin, ymax, ydiv, optmin, optmax, optstp, optwi);

      hist->GetYaxis()->SetRangeUser(ymin, ymax);
      if ( (ymax - optmax) < optwi/2 )
        hist->GetYaxis()->ChangeLabel( -1, -1, 0);
    }
  }

  //----------------------------------------------------------------------
  void PimpHist(TH1* hist, Style_t linestyle, Color_t linecolor, int linewidth, Style_t markerstyle, Color_t markercolor, double markersize){
    hist->SetLineColor(linecolor);
    hist->SetLineStyle(linestyle);
    hist->SetLineWidth(linewidth);
    hist->SetMarkerColor(markercolor);
    hist->SetMarkerStyle(markerstyle);
    hist->SetMarkerSize(markersize);
  }


  //----------------------------------------------------------------------
  void SplitCanvas(double ysplit, TPad*& p1, TPad*& p2)
  {
    if(!gPad) new TCanvas;
    p1 = new TPad("", "", 0, 0, 1, 1);
    p2 = new TPad("", "", 0, 0, 1, 1);

    p1->SetBottomMargin(ysplit);
    p2->SetTopMargin(1-ysplit);

    // Draw p1 second since it's often the more important one, that the user
    // would prefer to be able to interact with.
    for(TPad* p: {p2, p1}){
      p->SetFillStyle(0);
      p->Draw();
    }
  }


  //----------------------------------------------------------------------
  // split canvas in 4 for numu quartiles
  void SplitCanvasQuant(TCanvas *& canvas, TPad *& pad1, TPad *& pad2, TPad *& pad3, TPad *& pad4){

    canvas = new TCanvas(UniqueName().c_str(),"",960,800);

    pad3 = new TPad(UniqueName().c_str(),"pad3",0,0,1,1);
    pad3->Draw();
    pad3->SetTopMargin(0.5);
    pad3->SetBottomMargin(0.1);
    pad3->SetLeftMargin(0.1);
    pad3->SetRightMargin(0.49);
    pad3->SetFillStyle(0);
    canvas->cd();

    pad4 = new TPad(UniqueName().c_str(),"pad4",0,0,1,1);
    pad4->Draw();
    pad4->SetTopMargin(0.5);
    pad4->SetBottomMargin(0.1);
    pad4->SetLeftMargin(0.51);
    pad4->SetRightMargin(0.08);
    pad4->SetFillStyle(0);
    canvas->cd();

    pad1 = new TPad(UniqueName().c_str(),"pad1",0,0,1,1);// x1 y1 x2 y2
    pad1->Draw();
    pad1->SetTopMargin(0.1);
    pad1->SetBottomMargin(0.5);
    pad1->SetLeftMargin(0.1);
    pad1->SetRightMargin(0.49);
    pad1->SetFillStyle(0);
    canvas->cd();

    pad2 = new TPad(UniqueName().c_str(),"pad2",0,0,1,1);
    pad2->Draw();
    pad2->SetTopMargin(0.1);
    pad2->SetBottomMargin(0.5);
    pad2->SetLeftMargin(0.51);
    pad2->SetRightMargin(0.08);
    pad2->SetFillStyle(0);
    canvas->cd();

  }


  //----------------------------------------------------------------------
  void CenterTitles ( TH1 *  histo)
   {
     histo->GetXaxis()->CenterTitle();
     histo->GetYaxis()->CenterTitle();
     histo->GetZaxis()->CenterTitle();
   }


  //----------------------------------------------------------------------
  bool SortSystsName(const ISyst* systA, const ISyst* systB){
    return ( systA->ShortName() < systB->ShortName() );
  }


  TH1* PlotSystShifts(const SystShifts & shifts, bool sortName){

    auto systs = shifts.ActiveSysts();
    int nsysts = systs.size();

    // sort systematics by name
    if(sortName) std::sort(systs.begin(), systs.end(), SortSystsName);

    auto h = new TH1D ("sh",";; Pull (N#sigma)",
                       nsysts, -0.5, nsysts + 0.5);

    for (int systIdx = 0; systIdx  < nsysts; ++systIdx){
      double shiftThis = shifts.GetShift(systs[systIdx]);
      h->SetBinContent( systIdx + 1, shiftThis);

      TString tempName = systs[systIdx]->ShortName();
      if( tempName.Contains("Neutron") ) h->GetXaxis()->SetBinLabel( systIdx + 1, "NeutronSyst2018");
      else h->GetXaxis()->SetBinLabel( systIdx + 1, systs[systIdx]->ShortName().c_str());
    }

    h->SetMarkerStyle(kFullCircle);
    CenterTitles(h);
    h->SetTitleOffset(3.0);
    h->GetXaxis()->SetNdivisions(nsysts,kFALSE);
    h->GetXaxis()->SetLabelSize(0.065);
    h->GetXaxis()->LabelsOption("v");
    h->GetYaxis()->SetRangeUser(-1,1);
    auto c1 = new TCanvas ("c1","c1",700,350);
    c1->SetBottomMargin(0.5);
    h->Draw("lp hist");
    return h;
  }

  //----------------------------------------------------------------------
  TGraph* JoinGraphs(TGraph* a, TGraph* b, int fillcol)
  {
    TGraph* ret = (TGraph*)a->Clone(UniqueName().c_str());
    for(int i = 0; i < b->GetN(); ++i){
      double x, y;
      b->GetPoint(i, x, y);
      ret->SetPoint(ret->GetN(), x, y);
    }

    ret->SetFillColor(fillcol);

    return ret;
  }

  //----------------------------------------------------------------------
  TGraph* ExtendGraphToTop(TGraph* g, int col, double xmin, double xmax, double y)
  {
    g->SetPoint(g->GetN(), xmax, y);
    g->SetPoint(g->GetN(), xmin, y);

    g->SetFillColor(col);

    return g;
  }

  //----------------------------------------------------------------------
  TPaveText* DrawBeamLabel(bool isFHC)
  {

    TPaveText *pText = new TPaveText (0.1, 0.90, 0.16, 0.94,"NDC");
    pText->SetFillStyle(0);
    pText->SetLineColor(0);
    pText->SetLineWidth(0);
    pText->SetBorderSize(1);
    pText->SetTextColor(kGray+2);
    if (isFHC) pText->AddText("#nu-beam");
    else       pText->AddText("#bar{#nu}-beam"); 
    pText->SetTextSize(2/40.);  
    pText->SetTextAlign(11);
    pText->Draw();

    return pText;
  }

  //----------------------------------------------------------------------
  TPaveText* DrawQuantLabel(int quant)
  {
    TPaveText *ret = 0;
    // Draw nothing for q = 0
    if (quant == 0) {
      // Draw nothing for q = 0 but it's still valid
    }
    else if (quant == 1) {
      ret = new TPaveText(0.37, 0.75, 0.47, 0.90, "NDC");
      ret->SetBorderSize(0);
      ret->SetFillStyle(0);
      TText *text = ret->AddText(TString::Format("#splitline{#bf{Quartile %d}}{best resolution}", quant));
      text ->SetTextAlign(22);
      text ->SetTextSize(0.025);
      ret->Draw();
    }
    else if (quant == 2) {
      ret = new TPaveText(0.75, 0.75, 0.85, 0.90, "NDC");
      ret->SetBorderSize(0);
      ret->SetFillStyle(0);
      TText *text = ret->AddText(TString::Format("#bf{Quartile %d}", quant));
      text ->SetTextAlign(22);
      text ->SetTextSize(0.025);
      ret->Draw();
    }
    else if (quant == 3) {
      ret = new TPaveText(0.37, 0.35, 0.47, 0.50, "NDC");
      ret->SetBorderSize(0);
      ret->SetFillStyle(0);
      TText *text = ret->AddText(TString::Format("#bf{Quartile %d}", quant));
      text ->SetTextAlign(22);
      text ->SetTextSize(0.025);
      ret->Draw();
    }
    else if (quant == 4) {
      ret = new TPaveText(0.75, 0.35, 0.85, 0.50, "NDC");
      ret->SetBorderSize(0);
      ret->SetFillStyle(0);
      TText *text = ret->AddText(TString::Format("#splitline{#bf{Quartile %d}}{worst resolution}", quant));
      text ->SetTextAlign(22);
      text ->SetTextSize(0.025);
      ret->Draw();
    }
    else {
      std::cout << "\n***********************************************\n"
                << "\n  No valid quantile chosen for DrawQuantLabel  \n"
                << "\n***********************************************\n"
                << std::endl;
    }

    return ret;
  }

} // namespace