DrawUtils.h

Go to the documentation of this file.

#include "TH1.h"
#include "TH2.h"
#include "TMarker.h"
#include "TProfile.h"
#include "TProfile2D.h"
#include "TLine.h"
#include "TLegend.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TText.h"


typedef struct Plot Plot;
struct Plot {
  std::string  file;
  std::string  name;
  TH1*         h;
  std::string  type;
  std::string  legLabel;
  int          color;
  Plot(){
    h=0;
  }
  Plot(TH1* h_){
    h=h_;
    if(h)
      type=h->ClassName();
  }
  Plot(TH1* h_, std::string name_, 
       std::string legLabel_){
    h=h_;
    name=name_;
    legLabel=legLabel_;
    if(h)
      type=h->ClassName();
  }
  Plot(std::string file_, 
       std::string name_,
       std::string legLabel_){
    file=file_;
    name=name_;
    legLabel=legLabel_;
    TFile * f = TFile::Open(file.c_str());
    f->GetObject(name.c_str(),h);
    if(!h){ 
      std::cout << "did not find " << name
                << " in " << f->GetName() 
                << std::endl;
      std::abort();
    } else {
      h->SetDirectory(0);
      type=h->ClassName();
      if(h->GetEntries()==0){
        std::cout << "Zero entries in " << name
                  << " from " << f->GetName() 
                  << std::endl;
        //std::abort(); // occasionally okay to have zero entries if not ratio denom
      }
    }
    f->Close();
  }
};


class DrawUtils{

public:

  DrawUtils();
  ~DrawUtils();

  // Stacks plot_vec colored with "plc pmc" onto canv.
  //   If hDenom is specified, it is drawn first, grey and dotted.
  //   If ratio_title is provided, the ratios of 
  //     plot_vec over hDenom are drawn with the same colors.
  //   Instead of drawing a legend on canv, potentially obscuring info,
  //     draw a legend on canvleg for the user to place manually
  //   Plot is a class with the TH1*, label, etc..
  void DrawPlotStack( TCanvas *&canv, TCanvas *&canvleg,
                      const std::vector< Plot > &plot_vec, 
                      const Plot &plot_denom,
                      std::string ratio_title=""); // default empty
  void DrawPlotStack( TCanvas *&canv, TCanvas *&canvleg,
                      const std::vector< Plot > &plot_vec );
  
  void MakeRatio( TPad *&pad,
                  const std::vector<Plot> &plot_vec, 
                  const Plot &plot_denom, 
                  std::string ratio_name );
  void SetRatioErrors(TGraphErrors* &gRat, TH1* hNum, TH1* hDenom);

  void Draw2DPlot( TCanvas *&canv, TH2* h2 );

  void MakeLegend( TCanvas *&cleg, const std::vector< Plot > &plots, const Plot &denom);

  void GetBestYRange( float & ymin, float & ymax, 
                      std::vector< TProfile * > profiles,
                      float maxpad,
                      int rebin = 0, 
                      float outlierCut = 0.);
  void GetBestYRange( float & ymin, float & ymax, 
                      const std::vector< Plot > &plots,
                      float maxpad, int rebin = 0, float outlierCut = 0.);

  float GetLastYMin(){return fLastYMin;}
  float GetLastYMax(){return fLastYMax;}

private:

  int Cnum(){return fCnum; fCnum++;}
  int Pnum(){return fPnum; fPnum++;}

  int fCnum; //unique canvas number
  int fPnum; //unique plot nummber

  float fLastYMin;
  float fLastYMax;

};


//^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
DrawUtils::DrawUtils()
  : fCnum(0)
  , fPnum(0)
{
}
//^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
DrawUtils::~DrawUtils()
{
}


//^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void DrawUtils::DrawPlotStack( TCanvas *&c, TCanvas *&cleg,
                               const std::vector< Plot > &plot_vec )
{
  DrawPlotStack( c, cleg, plot_vec, Plot(0,"null","null plot"), "" );
}


//^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void DrawUtils::DrawPlotStack( TCanvas *&c, TCanvas *&cleg,
                               const std::vector< Plot > &plot_vec, 
                               const Plot &denom,
                               std::string ratio_title )
{
  // Make legend seperately, placing on plots automatically
  //  potentially covers up something important 
  MakeLegend(cleg,plot_vec, denom);
  c->Clear(); c->cd();
  gStyle->SetOptStat(0);

  // temporary until debugged for TH1F
  if(denom.type!="TProfile") return;

  std::string plotname;
  bool isProf;
  if(plot_vec.size()==0){
    std::cout << "DrawUtils::Draw Requires at least one plot" << std::endl;
    return;
  } else {
    plotname = plot_vec[0].name;
    isProf = plot_vec[0].type=="TProfile";
    if( denom.h &&
        plot_vec[0].type != denom.type ){
      std::cout << "Denom type: "<< plot_vec[0].type 
                <<" not the same as plot vector: "
                << denom.type << std::endl;
      return;
    }
  }


  // If ratio left blank, still draw denom_plot as
  //   grey in the main plot, but don't draw ratio
  bool drawRatio = (ratio_title!="");


  // Prepare pads for optional ratio plots
  TPad *p1 = new TPad("p1","p1", 0, 0, 1, 1); 
  TPad *p2 = new TPad("p2","p2", 0, 0, 1, 1); 
  p1->SetBottomMargin(0.45);
  p2->SetTopMargin(0.55);
  p1->SetFillStyle(0); p2->SetFillStyle(0);
  if(drawRatio) p1->cd();
  else c->cd();


  // AutoZooming on Y
  float ymin(0), ymax(0);
  int rebin(0);
  float outlierCut(0.);
  if( plotname.find("_time_") != std::string::npos ){
    rebin = 0; outlierCut = 10.; } 
  // TODO: Make GetBestYRange work with rebin, change line above to use it
  std::vector< Plot > all_plot_vec = plot_vec;
  all_plot_vec.push_back(denom);
  if( isProf ) GetBestYRange( ymin, ymax,
                              all_plot_vec,
                              0.1, rebin, outlierCut );
  //bool autozoom = false;
  bool autozoom = ( ymin != std::numeric_limits<float>::max() &&
                    ymax != std::numeric_limits<float>::min() );


  // Should Eventually remove name-specific tweaks like this
  if( plotname.find("_time_") != std::string::npos ){
    ymin = ymin - 0.2*(ymax - ymin);
    fLastYMin = ymin;
  }

  // 1d plots: Overlay all samples for this plot
  //  If "denom" sample is specified, draw it first
  //  Grey dotted line for denom sample (MC or control sample or...)
  bool firstdraw=true;
  if(denom.h){
    if(drawRatio) denom.h->GetXaxis()->SetLabelSize(0.);
    denom.h->SetLineWidth(1);
    denom.h->SetLineColor(kGray+2);
    denom.h->SetMarkerColor(kGray+2);
    if(autozoom){
      denom.h->SetMinimum(ymin);
      denom.h->SetMaximum(ymax);
    }
    denom.h->Draw( Form( "%s", isProf ? "" : "hist" ));
    firstdraw=false;
  }
  

  // Now the rest of the plots
  for( auto& plot: plot_vec ){
    if( !plot.h ){
      std::cout << "Plot " << plot.name 
                << " is null, skipping" << std::endl;
      return;
    }
    if( plot.h==denom.h ) continue; // already drew
    
    if( isProf ){
      plot.h->SetLineWidth(2);
      //if(rebin) plot.h->Rebin(rebin);
      if(autozoom){
        plot.h->SetMinimum(ymin);
        plot.h->SetMaximum(ymax);
      }
    } else {
      plot.h->Scale(1/plot.h->Integral());
      plot.h->GetYaxis()->SetTitle("Arb. Units");
    }
    plot.h->SetStats(0);

    //if(!isProf) PrintUnderOverFlow(h,plotname);

    if(firstdraw){
      if(drawRatio) plot.h->GetXaxis()->SetLabelSize(0.);
      plot.h->Draw(Form("%s plc pmc", isProf ? "" : "hist"));
      firstdraw=false;
    } else {
      plot.h->Draw(Form("%s same plc pmc", isProf ? "" : "hist"));
    }
    
  } // plot_vec loop

  
  if(drawRatio){
    MakeRatio( p2, plot_vec, denom, ratio_title );
    c->cd();
    p1->Draw();
    p2->Draw();
  }

  
} // DrawPlotStack


//^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void DrawUtils::Draw2DPlot( TCanvas *&canv,
                            TH2* h2 )
{
  std::string plotname;
  if(!h2->InheritsFrom(TH2::Class())){ 
    std::cout << "not a TH2: " << plotname << std::endl; return; }

  canv->cd();
  h2->SetStats(0);
  gStyle->SetPalette(kBird);
  h2->Draw("colz");
  //gROOT->SetStyle("novaStyle");
}




//^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void DrawUtils::MakeRatio( TPad *&pad, 
                           const std::vector<Plot> &plot_vec, 
                           const Plot &plot_denom, 
                           std::string ratio_title )
{
  // Build vector of ratios
  std::vector<TProfile*> ratio_vec;
  if(plot_denom.h->InheritsFrom(TProfile::Class()))
    for(auto& pl : plot_vec)
      ratio_vec.push_back((TProfile*)pl.h->Clone(Form("%s_%i",pl.h->GetName(),Pnum())));
  //if(plot_denom.h->InheritsFrom(TH1F::Class()))
  //  for(auto& pl : plot_vec)
  //    ratio_vec.push_back((TH1F*)pl.h->Clone(Form("%s_%i",pl.h->GetName(),Pnum())));
  if( plot_denom.type == plot_vec[0].type ){
    for(auto& rat : ratio_vec)
      rat->Divide(plot_denom.h);
  } else { 
    std::cout << "Denom different type, cannot take ratio" << std::endl;
    return;
  }

  pad->cd();
  TH1 *hAxes = (TH1*)plot_vec[0].h->Clone();
  hAxes->Draw("AXIS");
  hAxes->GetYaxis()->SetTitle(ratio_title.c_str());

  for( unsigned int i=0; i<ratio_vec.size(); i++ ){

    TGraphErrors *gRat = new TGraphErrors(ratio_vec[i]);
    gRat->SetLineColor  (plot_vec[i].h->GetLineColor());
    gRat->SetMarkerColor(plot_vec[i].h->GetLineColor());
    SetRatioErrors(gRat, plot_vec[i].h, plot_denom.h);

    if(i==0){
      gRat->Draw("P PLC PMC");
      double x[2] = {gRat->GetXaxis()->GetXmin(), gRat->GetXaxis()->GetXmax()};
      double y[2] = {1.0, 1.0};
      TGraph* l_Xaxis = new TGraph(2, x, y);
      l_Xaxis->SetLineColor(kGray);
      l_Xaxis->SetLineWidth(3);
      l_Xaxis->SetLineStyle(2);
      l_Xaxis->Draw("L SAME");
    }
    else {
      gRat->Draw("P PLC PMC");
    } 
  } // ratio_vec loop

  float ymin, ymax;
  GetBestYRange( ymin, ymax, ratio_vec, 0.1, 0, 10);
  hAxes->GetYaxis()->CenterTitle();
  hAxes->GetYaxis()->ChangeLabel();
  hAxes->GetYaxis()->SetRangeUser(ymin, ymax);

  return;
} // MakeRatio


//^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void DrawUtils::SetRatioErrors(TGraphErrors* &gRat, TH1* hNum, TH1* hDenom)
{
  for(int bin = 1; bin <= hNum->GetNbinsX(); bin++)
  {
    double num, denom, errnum, errdenom;
    num   = hNum->GetBinContent(bin);   errnum   = hNum->GetBinError(bin);
    denom = hDenom->GetBinContent(bin); errdenom = hDenom->GetBinError(bin);

    if(num!=0. && denom!=0.){
      gRat->SetPointError(bin-1, 0., (num/denom)*std::sqrt((errnum/num)*(errnum/num)+(errdenom/denom)*(errdenom/denom)));
    }
    else{
      gRat->SetPointError(bin-1, 0., 0.);
    }
  }

  return;
} // SetRatioErrors




//^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void DrawUtils::MakeLegend( TCanvas *&cleg, const std::vector< Plot > &plots, const Plot &denom)
{
  cleg->cd();
  TLegend *leg   = new TLegend(0, 0, 1, 1);
  if(plots.size()+1>2 && plots.size()+1<7) leg->SetNColumns(3);
  if(plots.size()+1>=7) leg->SetNColumns(5);

  leg->AddEntry(denom.h, denom.legLabel.c_str(), "L");
  for(auto &p : plots)
    leg->AddEntry(p.h, p.legLabel.c_str(), "L"); 

  leg->Draw();
  return;
} // MakeLegend



//^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void DrawUtils::GetBestYRange( float & ymin, float & ymax,
                               const std::vector< Plot > &plots,
                               float maxpad, int rebin, float outlierCut )
{
  if(plots[0].type!="TProfile"){
    std::cout << "GetBestYRange only works for profiles, skip" 
              << plots[0].name
              << std::endl;

    return;
  }
  std::vector< TProfile* > profiles;
  for( auto& pl : plots )
    profiles.push_back((TProfile*)pl.h);

  GetBestYRange( ymin, ymax, profiles, maxpad, rebin, outlierCut );
}



//^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void DrawUtils::GetBestYRange( float & ymin, float & ymax,
                               std::vector< TProfile* > profiles,
                               float maxpad,
                               int rebin, 
                               float outlierCut )
{
  ymin = std::numeric_limits<float>::max();
  ymax = std::numeric_limits<float>::min();
  TProfile * h;
  for( auto& p: profiles ){
    h = (TProfile *)p->Clone(Form("%s_%i",p->GetName(),Pnum()));
    //if(rebin) h->Rebin(rebin);

    for( int b = 1; b < h->GetXaxis()->GetNbins()-1; b++ ){
      // TProfile returns content=0 for bins with no entries.
      // Number of entries isn't publicly accessible, but the
      //   bin error is always 0 when there are no entries.
      if( h->GetBinError(b) == 0 ) continue;

      if( h->GetBinContent(b) < ymin ) ymin = h->GetBinContent(b);
      if( ymax < h->GetBinContent(b) ) ymax = h->GetBinContent(b);
    }
  }  


  if(outlierCut!=0.){    

    // Get RMS of all bins from all profiles combined by filling a TH1
    TH1F * binContDist = new TH1F(Form("%s_bin_contents_%i",
                                       profiles[0]->GetName(), Pnum()),"",
                                  1000,ymin,ymax);
    for( auto& p: profiles ){
      h = (TProfile *)p->Clone(Form("%s_%i",p->GetName(),Pnum()));
      //if(rebin) h->Rebin(rebin);

      for( int b = 1; b < h->GetXaxis()->GetNbins()-1; b++ ){
        if( h->GetBinError(b) == 0 ) continue;
        binContDist->Fill(h->GetBinContent(b));
      }
    }
    float combinedMean = binContDist->GetMean();
    float combinedRMS  = binContDist->GetRMS();
    //std::cout << "Combined Bins: mean = " << combinedMean 
    //          << ", rms = " << combinedRMS << std::endl;

    // Last loop cutting outliers using combined RMS
    ymin = std::numeric_limits<float>::max();
    ymax = std::numeric_limits<float>::min();    
    for( auto& p: profiles ){
      h = (TProfile *)p->Clone(Form("%s_%i",p->GetName(),Pnum()));
      //if(rebin) h->Rebin(rebin);
      for( int b = 1; b < h->GetXaxis()->GetNbins()-1; b++ ){
        if( h->GetBinError(b) == 0 ) continue;

        if( outlierCut*combinedRMS < std::abs(h->GetBinContent(b) - combinedMean) ) continue;

        if( h->GetBinContent(b) < ymin ) ymin = h->GetBinContent(b);
        if( ymax < h->GetBinContent(b) ) ymax = h->GetBinContent(b);
      }
    }  

  }

  if( ymin == std::numeric_limits<float>::max() ||
      ymax == std::numeric_limits<float>::min() ){
    std::cout << "Min/Max Y Failed in " << profiles[0]->GetName()
              << ": min = " << ymin << ", max = " << ymax 
              << "\n The plot may be empty."
              << std::endl;
    return;
  }

  float ymid = (ymax+ymin)/2;
  float fracRange = (ymax-ymin)/ymid;

  float pad = ( (maxpad - 0.1) / (0.01 - 0.2) ) * fracRange;
  if( pad > maxpad ) pad = maxpad; // with ceiling
  if( pad < 0.1 ) pad = 0.1; // and floor

  ymax = ymax + pad*ymid;
  ymin = ymin - pad*ymid;

  fLastYMax = ymax;
  fLastYMin = ymin;

  if(ymax < ymin){
    std::cout << "Min/Max Y Failed in " << profiles[0]->GetName()
              << ": min = " << ymin << ", max = " << ymax 
              << std::endl;
    return;
  }

} // GetBestYRange