Plotting.h

Go to the documentation of this file.

#include "CAFAna/XSec/ICrossSectionAnalysis.h"

#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TCanvas.h"
#include "TLegend.h"
namespace ana {

  /// General plotting function that plots the vector of hists
  /// to a multi-page pdf, where the first page can be an overlay of
  /// all of the hists
  void PlotAll(TString outfile,
               TString title,
               std::vector<TH1*> hists,
               std::vector<TString> labels,
               TString draw_option,
               bool overlay,
               bool logx = false)
  {
    TLegend * leg = new TLegend();
    TCanvas * c = new TCanvas();
    c->Print(outfile+"["); //opens pdf for printing
    c->cd();
    double maxy = -1;
    for(int i = 0; i < (int) hists.size(); i++) {
      hists[i]->GetXaxis()->CenterTitle();
      hists[i]->GetYaxis()->CenterTitle();
      if(hists[i]->GetMaximum() > maxy) maxy = hists[i]->GetMaximum();
    }
    
    // leave some headroom
    maxy = 1.2 * maxy;
    
    if(overlay) {
      hists[0]->SetTitle(TString::Format("%s: All",
                                         title.Data()));
      hists[0]->GetYaxis()->SetRangeUser(0, maxy);
      hists[0]->Draw(draw_option);
      leg->AddEntry(hists[0], labels[0]);
      for(int i = 1; i < (int) hists.size(); i++) {
        hists[i]->Draw(draw_option + "same");
        leg->AddEntry(hists[i], labels[i]);
      }
      leg->Draw();
      if(logx) c->SetLogx();
      c->Print(outfile);
    }
    for(int i = 0; i < (int) hists.size(); i++) {
      hists[i]->SetTitle(TString::Format("%s: %s",
                                         title.Data(),
                                         labels[i].Data()));
      hists[i]->GetYaxis()->SetRangeUser(0, maxy);
      hists[i]->Draw(draw_option);
      if(logx) c->SetLogx();
      c->Print(outfile);
    }
    c->Print(outfile+"]");
  }

  void PlotAllSelectionDecomposition(TString outfile,
                                     ICrossSectionAnalysis * xsec_ana,
                                     std::vector<XSecSysts::Syst_t> systs,
                                     std::vector<TString> systs_labels,
                                     std::vector<TString> bkgd_labels,
                                     std::vector<Int_t> bkgd_colors)
  {
    TCanvas * c = new TCanvas();
    c->Print(outfile+"["); //opens pdf for printing
    c->cd();
    for(int isyst = 0; isyst < (int) systs.size(); isyst++) {
      int ncomponents = xsec_ana->GetCrossSectionSpectra(systs[isyst])->GetSignalEstimator()->GetNBkgdComponents();
      std::vector<TH1*> bkgds(ncomponents);
      TLegend * leg = new TLegend();
      
      TH1 * selected = xsec_ana->GetCrossSectionSpectra(systs[isyst])->GetSelectedAllRecoSpace()->ToTH1(xsec_ana->GetCrossSectionSpectra(systs[isyst])->GetSelectedAllRecoSpace()->POT());
      selected->SetTitle("Selection decomposition: " + systs_labels[isyst]);
      selected->Draw("hist");
      leg->AddEntry(selected, "Selected", "l");
      for(int ibkgd = 0; ibkgd < ncomponents; ibkgd++) {
        bkgds[ibkgd] = (TH1*) xsec_ana->GetCrossSectionSpectra(systs[isyst])->GetSignalEstimator()->BackgroundComponent(ibkgd)->ToTH1(xsec_ana->GetCrossSectionSpectra(systs[isyst])->GetSignalEstimator()->BackgroundComponent(ibkgd)->POT())->Clone();
        bkgds[ibkgd]->SetLineColor(bkgd_colors[ibkgd]);
        bkgds[ibkgd]->Draw("hist same");
        leg->AddEntry(bkgds[ibkgd], bkgd_labels[ibkgd], "l");   
      }
      leg->Draw();
      c->Print(outfile);
    }    
    c->Print(outfile+"]");
  }

  /// Plot all efficiencies
  void PlotAllEfficiency(TString outfile,
                         ICrossSectionAnalysis * xsec_ana,
                         std::vector<XSecSysts::Syst_t> systs,
                         std::vector<TString> labels,
                         std::vector<Int_t> colors)
  {
    std::vector<TH1*> eff;
    for(int isyst = 0; isyst < (int) systs.size(); isyst++) {
      eff.push_back((TH1*) xsec_ana->Efficiency(systs[isyst])->Clone());
      eff.back()->SetLineColor(colors[isyst]);
      eff.back()->GetYaxis()->SetTitle("#epsilon");
      eff.back()->GetYaxis()->CenterTitle();
    }
    PlotAll(outfile,
            "Efficiency",
            eff,
            labels,
            "hist",
            true);
  }

  /// Plot all flux spectra
  void PlotAllFlux(TString outfile,
                   ICrossSectionAnalysis * xsec_ana,
                   std::vector<XSecSysts::Syst_t> systs,
                   std::vector<TString> labels,
                   std::vector<Int_t> colors)
  {
    std::vector<TH1*> flux;
    for(int isyst = 0; isyst < (int) systs.size(); isyst++) {
      flux.push_back(xsec_ana->Flux(systs[isyst]));
      flux.back()->SetLineColor(colors[isyst]);
    }
    PlotAll(outfile,
            "Flux",
            flux,
            labels,
            "hist",
            true);
  }

  /// Plot all reco to true matrices
  void PlotAllRecoTrue(TString outfile,
                       ICrossSectionAnalysis * xsec_ana,
                       std::vector<XSecSysts::Syst_t> systs,
                       std::vector<TString> labels)

  {
    std::vector<TH1*> recotrue;
    for(int isyst = 0; isyst< (int) systs.size(); isyst++) {
      double POT = xsec_ana->GetCrossSectionSpectra(systs[isyst])->GetRecoTrue()->POT();
      recotrue.push_back(xsec_ana->GetCrossSectionSpectra(systs[isyst])->GetRecoTrue()->ToTH2(POT));
    }

    PlotAll(outfile,
            "Unfolding Matrix",
            recotrue,
            labels,
            "colz",
            false);    
  }

  /// Plot all of the signal estimates in a ICrossSectionAnalysis object
  void PlotAllSignalEstimates(TString outfile,
                              ICrossSectionAnalysis * xsec_ana,
                              std::vector<XSecSysts::Syst_t> systs,
                              std::vector<TString> labels,
                              std::vector<Int_t> colors)
  {
    std::vector<TH1*> signalest;
    for(int isyst = 0; isyst < (int) systs.size(); isyst++) {
      double POT = xsec_ana->SignalEstimate(systs[isyst])->POT();
      signalest.push_back(xsec_ana->SignalEstimate(systs[isyst])->ToTH1(POT));
      signalest.back()->SetLineColor(colors[isyst]);
      signalest.back()->SetTitle(TString::Format("Signal Esimate: %s",
                                                 labels[isyst].Data()));
    }
      
    PlotAll(outfile,
            "Signal Estimate",
            signalest,
            labels,
            "hist",
            true);          
  }

  /// Plot all unfolded signal estimates in an ICrossSectionAnalysis object
  void PlotAllUnfoldedSignalEstimates(TString outfile,
                                      ICrossSectionAnalysis * xsec_ana,
                                      std::vector<XSecSysts::Syst_t> systs,
                                      std::vector<TString> labels,
                                      std::vector<Int_t> colors)
  {    
    std::vector<TH1*> signal;
    for(int isyst = 0; isyst < (int) systs.size(); isyst++) {
      signal.push_back((TH1*) xsec_ana->UnfoldedSignal(systs[isyst])->Clone());
      signal.back()->SetLineColor(colors[isyst]);
    }
      
    PlotAll(outfile,
            "Unfolded Signal",
            signal,
            labels,
            "hist",
            true);          
  }

  /// Plot all purity distributions in an ICrossSectionAnalysis object
  void PlotAllPurity(TString outfile,
                     ICrossSectionAnalysis * xsec_ana,
                     std::vector<XSecSysts::Syst_t> systs,
                     std::vector<TString> labels,
                     std::vector<Int_t> colors)
  {
    std::vector<TH1*> purity;
    for(int isyst = 0; isyst < (int) systs.size(); isyst++) {
      purity.push_back((TH1*) xsec_ana->Purity(systs[isyst])->Clone());
      purity.back()->SetLineColor(colors[isyst]);
    }

    PlotAll(outfile,
            "Purity",
            purity, 
            labels,
            "hist",
            true);
  }

  /// Plot cross section over Enu for each systematic shift
  void PlotAllCrossSectionOverENu(TString outfile,
                                  ICrossSectionAnalysis * xsec_ana,
                                  std::vector<XSecSysts::Syst_t> systs,
                                  std::vector<TString> syst_labels,
                                  std::vector<Int_t> colors,
                                  bool logx = true)
  {
    std::vector<TH1*> xsec_over_e;
    for(int isyst = 0; isyst < (int) systs.size(); isyst++) {
      xsec_over_e.push_back((TH1*) xsec_ana->GetCrossSection(systs[isyst])->Clone());
      xsec_over_e.back()->SetLineColor(colors[isyst]);
      xsec_over_e.back()->GetXaxis()->SetTitle("Neutrino Energy (GeV)");
      xsec_over_e.back()->GetYaxis()->SetTitle("#sigma / E_{#nu} (10^{-38} cm^{2} / GeV)");
      xsec_over_e.back()->GetXaxis()->CenterTitle();
      xsec_over_e.back()->GetYaxis()->CenterTitle();
      for(int ibin = 1; ibin <= xsec_over_e.back()->GetNbinsX(); ibin++) {
        double e = xsec_over_e.back()->GetBinCenter(ibin);
        double xsec = xsec_over_e.back()->GetBinContent(ibin);
        xsec_over_e.back()->SetBinContent(ibin, xsec / e);
      }
      xsec_over_e.back()->Scale(1/1e-38);
    }
    PlotAll(outfile,
            "#sigma / E_{#nu}",
            xsec_over_e,
            syst_labels,
            "hist",
            true,
            logx);
  }
  
  /// Plot cross section result of each systematic shift
  void PlotAllCrossSectionResult(TString outfile,
                                 TString xtitle,
                                 TString ytitle,
                                 ICrossSectionAnalysis * xsec_ana,
                                 std::vector<XSecSysts::Syst_t> systs,
                                 std::vector<TString> syst_labels,
                                 std::vector<Int_t> colors,
                                 double scale = 1e-38)
  {
    std::vector<TH1*> result;
    for(int isyst = 0; isyst < (int) systs.size(); isyst++) {
      result.push_back((TH1*) xsec_ana->GetCrossSection(systs[isyst])->Clone());
      result.back()->SetLineColor(colors[isyst]);
      if(scale > 0) result.back()->Scale(1 / scale);
      result.back()->GetXaxis()->SetTitle(xtitle.Data());
      result.back()->GetYaxis()->SetTitle(ytitle.Data());
    }
    PlotAll(outfile,
            "Cross Section Result",
            result,
            syst_labels,
            "hist",
            true);          
  }

  TGraphAsymmErrors * GetAsymmErrorBars(TH1D * hnom, TH1D * hup, TH1D * hdown)
  {
    TGraphAsymmErrors * g = new TGraphAsymmErrors();
    for(int ibin = 1; ibin <= hnom->GetNbinsX(); ibin++) {
      double y = hnom->GetBinContent(ibin);
      double x = hnom->GetBinCenter(ibin);
      double w = hnom->GetBinWidth(ibin);
      double err_up   = hup->GetBinContent(ibin);
      double err_down = hdown->GetBinContent(ibin);
      g->SetPoint(ibin, x, y);
      g->SetPointError(ibin, w/2, w/2, err_up, err_down);      
    }
    return g;
  }


}