Plot.cxx

Go to the documentation of this file.

// $Id: Plot.cxx,v 1.6 2007/09/13 21:23:56 rustem Exp $

// C++
#include <algorithm>
#include <cassert>
#include <cmath>
#include <fstream>
#include <iostream>
#include <iomanip>
#include <sstream>

// ROOT
#include "TArrow.h"
#include "TBox.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TGraphAsymmErrors.h"
#include "TH1.h"
#include "TF1.h"
#include "TLine.h"
#include "TPaveStats.h"
#include "TStyle.h"

// Local
#include "PhysicsNtuple/DataBlock.h"
#include "PhysicsNtuple/Hist1d.h"
#include "Plot.h"

using namespace std;

//-----------------------------------------------------------------------------------------
double Plot::GetChi2(const Hist &lhs, const Hist &rhs, const string &option, ostream &o)
{
   TH1 *hL = lhs.GetHist();
   TH1 *hR = rhs.GetHist();

   if(!hL || !hR)
   {
      cerr << "GetChi2() - null histogram pointers" << endl;
      return 0.0;
   }
   if(hL -> GetNbinsX() != hR -> GetNbinsX())
   {
      cerr << "GetChi2() - histograms have unequal number of bins" << endl;
      return 0.0;
   }

   unsigned int ndof = 0;
   double chi2 = 0.0;

   int lbin = 1, hbin = hL -> GetNbinsX();

   if(option.find("U") != string::npos || option.find("u") != string::npos)
   {
      --lbin;
   }
   if(option.find("O") != string::npos || option.find("o") != string::npos)
   {
      ++hbin;
   }
   if(!(lbin < hbin))
   {
      cerr << "GetChi2() - invalid range of histogram bins" << endl;
      return 0.0;
   }

   const map<int, pair<Band, Band> >& bandL = lhs.GetErrorBand();
   const map<int, pair<Band, Band> >& bandR = rhs.GetErrorBand();

   for(int ibin = lbin; ibin <= hbin; ++ibin)
   {
      const double valueL = hL -> GetBinContent(ibin);
      const double valueR = hR -> GetBinContent(ibin);
    
      if(!(valueL > 0.0) && !(valueR > 0.0))
      {
         continue;
      }
    
      const double nom = (valueL - valueR)*(valueL - valueR);
      double den = 0.0;

      const map<int, pair<Band, Band> >::const_iterator lit = bandL.find(ibin);
      if(lit != bandL.end())
      {
         const Band &band = (lit -> second).second;
         assert(band.Bin() == ibin && "Mismatched bin number");

         if(band.Low() > band.High())
         {
            den += band.Low();
         }
         else
         {
            den += band.High();
         }
      }
      else
      {
         den += valueL;
      }
      
      const map<int, pair<Band, Band> >::const_iterator rit = bandR.find(ibin);
      if(rit != bandR.end())
      {
         const Band &band = (rit -> second).second;
         assert(band.Bin() == ibin && "Mismatched bin number");

         if(band.Low() > band.High())
         {
            den += band.Low();
         }
         else
         {
            den += band.High();
         }
      }
      else
      {
         den += valueR;
      }

      if(!(den > 0.0))
      {
         cerr << "Plot::GetChi2 - bin error is zero" << endl;
         continue;
      }

      chi2 += nom/den;
      ++ndof;
   }

   if(!(ndof > 0))
   {
      return 0.0;
   }

   if(option.find("P") != string::npos || option.find("V") != string::npos)
   {
      o << "chi2/ndof = " << chi2 << "/" << ndof << " = " << chi2/ndof << endl;
   }

   if(option.find("chi2") != string::npos)
   {
      return chi2;
   }
   else if(option.find("ndof") != string::npos)
   {
      return ndof;
   }

   return chi2/ndof;
}

//-----------------------------------------------------------------------------------------
Plot::Style::Style()
   :lwidth(2),
    lstyle(1),
    lcolor(1),
    mstyle(-1),
    msize(0),
    mcolor(0),
    stat(1),
    rebin(0),
    xmin(+1.0),
    xmax(-1.0),
    statx(0.75),
    staty(0.5),
    statw(0.25),
    stath(0.23),
    option(),
    xtitle(),
    ytitle()
{
}

//-----------------------------------------------------------------------------------------
Plot::Style::~Style()
{
}

//-----------------------------------------------------------------------------------------
void Plot::Style::Print(ostream &o) const
{
   o << "lwidth = " << lwidth << endl;
   o << "lstyle = " << lstyle << endl;
   o << "lcolor = " << lcolor << endl;
   o << "mstyle = " << mstyle << endl;
   o << "msize = " << msize  << endl;
   o << "mcolor = " << mcolor << endl;
   o << "stat = " << stat   << endl;
   o << "rebin = " << rebin  << endl;
   o << "xmin = " << xmin << endl;
   o << "xmax = " << xmax << endl;
   o << "statx = " << statx << endl;
   o << "staty = " << staty << endl;
   o << "statw = " <<  statw << endl;
   o << "stath = " << stath << endl;
   o << "option = " << option << endl;
   o << "xtitle = " << xtitle << endl;
   o << "ytitle = " << ytitle << endl;
}

//-----------------------------------------------------------------------------------------
bool Plot::Style::Read(string path)
{
   std::ifstream infile(path.c_str());
   if(!infile || !infile.is_open())
   {
      cerr << "Style::Read - failed to open file " << path << endl;
      return false;
   }

   
   while(!infile.eof())
   {
      std::string line;
      std::getline(infile, line);

      const string::size_type epos = line.find("=");
      
      if(line.empty() || line.find("#") != string::npos || epos == string::npos)
      {
         continue;
      }

      Set(line.substr(0, epos), line.substr(epos + 1));
   }

   return true;
}

//-----------------------------------------------------------------------------------------
bool Plot::Style::Set(string name, string value)
{
   //cout << "(name, value) = (" << name << ", " << value << ")" << endl;

   if(name.find("lwidth") != string::npos) lwidth = atoi(value.c_str());
   if(name.find("lstyle") != string::npos) lstyle = atoi(value.c_str());
   if(name.find("lcolor") != string::npos) lcolor = atoi(value.c_str());
   if(name.find("mstyle") != string::npos) mstyle = atoi(value.c_str());
   if(name.find("msize" ) != string::npos) msize  = atoi(value.c_str());
   if(name.find("mcolor") != string::npos) mcolor = atoi(value.c_str());
   if(name.find("stat"  ) != string::npos) stat   = atoi(value.c_str());
   if(name.find("rebin" ) != string::npos) rebin  = atoi(value.c_str());
   if(name.find("xmin"  ) != string::npos) xmin   = atof(value.c_str());
   if(name.find("xmax"  ) != string::npos) xmax   = atof(value.c_str());
   if(name.find("statx" ) != string::npos) statx  = atof(value.c_str());
   if(name.find("staty" ) != string::npos) staty  = atof(value.c_str());
   if(name.find("statw" ) != string::npos) statw  = atof(value.c_str());
   if(name.find("stath" ) != string::npos) stath  = atof(value.c_str());
   if(name.find("option") != string::npos) option = value;
   if(name.find("xtitle") != string::npos) xtitle = value;
   if(name.find("ytitle") != string::npos) ytitle = value;

   return true;
}

//-----------------------------------------------------------------------------------------
double Plot::GetChi2(TH1 *h1, TH1 *h2, const string &option, const string &tag)
{
   if(!h1 || !h2)
   {
      cerr << "GetChi2() - invalid histogram pointers" << endl;
      return 0.0;
   }
   if(h1 -> GetNbinsX() != h2 -> GetNbinsX())
   {
      cerr << "GetChi2() - histograms have unequal number of bins" << endl;
      return 0.0;
   }

   double ndof = 0.0;
   double chi2 = 0.0;

   int lbin = 1, hbin = h1 -> GetNbinsX();

   if(option.find("U") != string::npos || option.find("u") != string::npos)
   {
      --lbin;
   }
   if(option.find("O") != string::npos || option.find("o") != string::npos)
   {
      ++hbin;
   }
   if(!(lbin < hbin))
   {
      cerr << "GetChi2() - invalid range of histogram bins" << endl;
      return 0.0;
   }

   for(int ibin = lbin; ibin <= hbin; ++ibin)
   {
      const double value1 = h1 -> GetBinContent(ibin);
      const double value2 = h2 -> GetBinContent(ibin);
    
      if(value1 > 0.0 && value2 > 0.0)
      {
         chi2 += (value1 - value2) * (value1 - value2)/value1;
         ndof += 1.0;
      }
   }

   if(!(ndof > 0.0))
   {
      return 0.0;
   }

   if(option.find("P") != string::npos || option.find("V") != string::npos || !tag.empty())
   {
      if(tag.empty())
      {
         cout << "chi2/ndof = " << chi2 << "/" << ndof << " = " << chi2/ndof << endl;
      }
      else
      {
         cout << tag << ": chi2/ndof = " << chi2 << "/" << ndof << " = " << chi2/ndof << endl;
      }
   }

   return chi2/ndof;
}

//-----------------------------------------------------------------------------------------
const string Plot::Format(const double v, const int w, const int p)
{
   stringstream ss;
   ss << std::fixed << std::setw(w) << std::setprecision(p) << v;
   return ss.str();
}

//-----------------------------------------------------------------------------------------
Plot::Hist::Hist(string fpath)
   :fFilePath(),
    fHistPath(),
    fHistPathL(),
    fHistPathH(),
    fHist(0),
    fHistL(0),
    fHistH(0),
    fProtons(0.0),
    fErrorColor(-1),
    fErrorStyle(1001),
    fAGraph(0),
    fError(),
    fArray()
{
   Set(fpath);
}

//-----------------------------------------------------------------------------------------
Plot::Hist::~Hist()
{
   if(fHist)  delete fHist;
   if(fHistL) delete fHistL;
   if(fHistH) delete fHistH;   

   ClearError();
}

//-----------------------------------------------------------------------------------------
bool Plot::Hist::Set(const string fpath)
{
   if(fpath.size() < 6)
   {
      return false;
   }

   TFile *file = new TFile(fpath.c_str(), "READ");
   if(!file -> IsOpen())
   {
      cerr << "Failed to open ROOT file " << fpath << endl;
      return false;
   }

   fFilePath = fpath;
   
   return true;
}

//-----------------------------------------------------------------------------------------
void Plot::Hist::SetProtons(const double protons)
{
   if(protons > 0.0)
   {
      fProtons = protons;
   }
}

//---------------------------------------------------------------------------------------
TH1* Plot::Hist::Get(const string hpath)
{   
   pair<TH1 *, double> result = Read(hpath);

   if(!result.first)
   {
      return 0;
   }
   
   fHistPath = hpath;

   fHist = result.first;

   if(result.second > 0.0)
   {
      fProtons = result.second;
   }

   return fHist;
}

//---------------------------------------------------------------------------------------
TGraphAsymmErrors* Plot::Hist::SetLH(const string lpath, const string hpath)
{   
   pair<TH1 *, double> lresult = Read(lpath);
   pair<TH1 *, double> hresult = Read(hpath);

   if(!lresult.first || !hresult.first)
   {
      return 0;
   }

   if(fHistL) delete fHistL;
   if(fHistH) delete fHistH;

   fHistPathL = lpath;
   fHistPathH = hpath;

   fHistL = lresult.first;
   fHistH = hresult.first;

   return ComputeError();
}

//---------------------------------------------------------------------------------------
TGraphAsymmErrors* Plot::Hist::ComputeError()
{   
   TH1 *hC = fHist;
   TH1 *hL = fHistL;
   TH1 *hH = fHistH;

   if(!hC || !hL || !hH)
   {
      return 0;
   }

   ClearError();

   const int nbin = hC -> GetNbinsX();

   if(hL -> GetNbinsX() != nbin || hH -> GetNbinsX() != nbin)
   {
      cerr << "Error::Compute - mismatched number of bins" << endl;
      return 0;
   }

   int npoint = 0;
   for(int ibin = 1; ibin <= nbin; ++ibin)
   {
      if(hC -> GetBinContent(ibin) > 0)
      {
         ++npoint;
      }
   }

   double* x   = Create(npoint);
   double* y   = Create(npoint);
   double* exl = Create(npoint);
   double* exh = Create(npoint);
   double* eyl = Create(npoint); 
   double* eyh = Create(npoint);

   npoint = 0;
   for(int ibin = 1; ibin <= nbin; ++ibin)
   {      
      if(!(hC -> GetBinContent(ibin) > 0))
      {
         continue;
      }

      const double xc = hC -> GetXaxis() -> GetBinCenter(ibin);
      const double xl = hC -> GetXaxis() -> GetBinLowEdge(ibin);
      const double xh = hC -> GetXaxis() -> GetBinUpEdge(ibin);

      const double yc = hC -> GetBinContent(ibin);
      const double yl = hL -> GetBinContent(ibin);
      const double yh = hH -> GetBinContent(ibin);

      if(yl > yh)
      {
         cerr << "Error::Compute -  low value " << yl << " larger than high value " << yh<< endl;
         continue;
      }

      if(yc < yl || yh < yc)
      {
         cerr << "Error::Compute -  low and/or high error have wrong limit" << yh << endl;
         continue;
      }

      const double error_stat = std::sqrt(yc);
      const double error_low  = error_stat + yc - yl;
      const double error_high = error_stat + yh - yc;

      //cout << "yc - yl = " << yc << " - " << yl << " = " << yc - yl << endl;
      //cout << "yh - yc = " << yh << " - " << yc << " = " << yh - yc << endl;

      x[npoint] = xc;
      y[npoint] = yc;

      exl[npoint] = xc - xl;
      exh[npoint] = xh - xc;

      eyl[npoint] = error_low;
      eyh[npoint] = error_high;

      ++npoint;

      const Band xb(ibin, xc, xl, xh);
      const Band yb(ibin, yc, yl - error_stat, yh + error_stat);
      
      const pair<Band, Band> np(xb, yb);

      if(!fError.insert(map<int, pair<Band, Band> >::value_type(ibin, np)).second)
      {
         cerr << "Error::Compute - bin " << ibin << " already exists" << endl;
      }

      TBox *box = new TBox(xb.Low(), yb.Low(), xb.High(), yb.High());

      if(!fBox.insert(map<int, TBox *>::value_type(ibin, box)).second)
      {
         cerr << "Error::Compute - bin " << ibin << " already exists" << endl;
      }
      
      
      if(error_low > error_high)
      {
         hC -> SetBinError(ibin, error_low);
      }
      else
      {
         hC -> SetBinError(ibin, error_high);
      }
   }

   fAGraph = new TGraphAsymmErrors(npoint, x, y, exl, exh, eyl, eyh);

   return fAGraph;
}

//---------------------------------------------------------------------------------------
TH1* Plot::Hist::operator->() const
{
   return fHist;
}

//---------------------------------------------------------------------------------------
double Plot::Hist::GetEntries(string option) const
{
   TH1 *h = fHist;

   if(option.find("l") != string::npos)
   {
      h = fHistL;
   }
   else if(option.find("h") != string::npos)
   {
      h = fHistH;
   }

   if(!h)
   {
      return 0.0;
   }

   long double area = 0.0;

   int lbin = 1, hbin = h -> GetNbinsX();

   if(option.find("u") != string::npos)
   {
      --lbin;
   }
   if(option.find("o") != string::npos)
   {
      ++hbin;
   }

   for(int ibin = lbin; ibin <= hbin; ++ibin)
   {
      const long double value = h -> GetBinContent(ibin);

      if(value > 0.0)
      {
         area += value;
      }
   }

   return static_cast<double>(area);
}

//---------------------------------------------------------------------------------------
double Plot::Hist::GetProtons() const
{
   return fProtons;
}

//---------------------------------------------------------------------------------------
const map<int, pair<Plot::Band, Plot::Band> >& Plot::Hist::GetErrorBand() const
{
   return fError;
}

//---------------------------------------------------------------------------------------
void Plot::Hist::Rebin(const int rebin)
{
   if(!fHist || rebin < 1)
   {
      return;
   }

   fHist -> Rebin(rebin);
   
   if(fHistL && fHistH)
   {
      fHistL -> Rebin(rebin);
      fHistH -> Rebin(rebin);
      
      ComputeError();
   }   
}

//---------------------------------------------------------------------------------------
TH1* Plot::Hist::Scale(const Hist &other, string option)
{
   if(!fHist)
   {
      return 0;
   }

   const double area = GetEntries("uo");

   double scale = 1.0;

   if(option.find("area") != string::npos)
   {
      if(area > 0.0)
      {
         scale = 1.0/area;
      }
   }
   else if(option.find("ent") != string::npos)
   {
      const double area_other = other.GetEntries("");
      if(area > 0.0 && area_other > 0.0)
      {
         scale = area_other/area;
      }
   }
   else if(option.find("pot") != string::npos)
   {
      if(fProtons > 0.0 && other.GetProtons() > 0.0)
      {
         scale = other.GetProtons()/fProtons;
      }
   }

   const double entries  = GetEntries("uo");
   const double entriesl = GetEntries("l uo");
   const double entriesh = GetEntries("h uo");

   fHist -> Scale(scale);
   
   if(fHistL) fHistL -> Scale(scale);
   if(fHistH) fHistH -> Scale(scale);

   if(option.find("area") == string::npos)
   {
      fHist -> SetEntries(entries * scale);
      
      if(fHistL) fHistL -> SetEntries(entriesl * scale);
      if(fHistH) fHistH -> SetEntries(entriesh * scale);
   }

   ComputeError();

   return fHist;
}

//---------------------------------------------------------------------------------------
TCanvas* Plot::Hist::Draw(string option, string cname)
{
   if(!fHist)
   {
      return 0;
   }
   
   TH1 *h = fHist;

   TCanvas *canvas = 0;

   if(!cname.empty())
   {
      canvas = new TCanvas(cname.c_str(), cname.c_str(), 0, 0, 650, 400);
      canvas -> cd();

      if(option.find("grid") != string::npos)
      {
         canvas -> SetGrid();
      }
      if(cname.find("log") != string::npos)
      {
         canvas -> SetLogy();
      }
      canvas -> Draw();
   }

   h -> Draw(option.c_str());

   if(option.find("E") != string::npos)
   {
      TH1 *hc = dynamic_cast<TH1 *> (fHist -> Clone("hc"));
      hc -> SetDirectory(0);
      hc -> Draw("same HIST");
   }

   if(fErrorColor > 0 && fErrorStyle > 0)
   {
      const int bmin = fHist -> GetXaxis() -> GetFirst();
      const int bmax = fHist -> GetXaxis() -> GetLast();

      for(map<int, TBox *>::const_iterator it = fBox.begin(); it != fBox.end(); ++it)
      {
         const int ibin = it -> first;
         
         if(ibin < bmin || ibin > bmax)
         {
            continue;
         }

         TBox *box = it -> second;
         if(!box)
         {
            continue;
         }
         
         box -> SetFillColor(fErrorColor);
         box -> SetFillStyle(fErrorStyle);       
         box -> Draw();
      }
   }

   gPad -> Modified();
   gPad -> Update();

   return canvas;
}

//---------------------------------------------------------------------------------------
void Plot::Hist::SetErrorStyle(const int color, const int style)
{
   fErrorColor = color;
   fErrorStyle = style;
}

//-----------------------------------------------------------------------------------------
void Plot::Hist::Move(float x, float y, float w, float h, string option)
{
   if(!fHist || !gPad)
   {
      return;
   }

   gPad -> Modified();
   gPad -> Update();

   TPaveStats* stat = dynamic_cast<TPaveStats *> (fHist -> FindObject("stats"));
   if(!stat)
   {
      if(option.find("V") != string::npos)
      {
         cerr << "Hist::Move - no statistics box is present" << endl;
      }
      return;
   }

   int color = fHist -> GetLineColor();
   if(option.find("M") != string::npos)
   {
      color = fHist -> GetMarkerColor();
   }

   stat -> SetStatFormat("6.3f");
   stat -> SetTextColor(color);
   stat -> SetLineColor(color);
   //stat -> SetOptStat(1111110);

   if(option.find("V") != string::npos)
   {
      cout << "Hist::Move - left lower corner = (" << x << ", " << y << ") and " 
           << "right upper corner = (" << x + w << ", " << y + h << ")" << endl;      
   }

   stat -> SetX1NDC(x);
   stat -> SetY1NDC(y);
   stat -> SetX2NDC(x + w);
   stat -> SetY2NDC(y + h);

   gPad -> Modified();
   gPad -> Update();
}


//-----------------------------------------------------------------------------------------
TH1* Plot::Hist::GetHist(string option) const
{
   TH1 *h = fHist;

   if(option.find("l") != string::npos)
   {
      h = fHistL;
   }
   else if(option.find("h") != string::npos)
   {
      h = fHistH;
   }

   return h;   
}

//-----------------------------------------------------------------------------------------
TH1* Plot::Hist::GetRatio(const Hist &other, string option,
                          const double min_, const double max_, const double factor) const
{
   if(!fHist || !other.GetHist())
   {
      cerr << "Error::GetRatio - invalid histogram pointer(s)" << endl;
      return 0;
   }

   TH1 *hC = fHist;
   TH1 *hL = fHistL;
   TH1 *hH = fHistH;
   TH1 *hO = other.GetHist();
   
   const int nbin = hC -> GetNbinsX();
   const double nmin = 1000;      
      
   if(hC -> GetNbinsX() != nbin)
   {
      cerr << "Error::GetRatio - mismatched number of bins" << endl;
      return 0;
   }
   
   TH1 *hR = CreateRatio(hC, hO, option, nmin);
   if(!hR)
   {
      cerr << "Error::GetRatio - failed to create ratio histogram" << endl;
      return 0;      
   }
   else if(option.find("merge") != string::npos)
   {
      if(hL && hH)
      {
         Anp::Hist1d<double> histL(*hR);
         Anp::Hist1d<double> histH(*hR);         

         histL.Reset();
         histH.Reset();
         
         histL.Fill(*hL);
         histH.Fill(*hH);

         hL = Anp::CreateTH1<double>(histL, "hL");
         hH = Anp::CreateTH1<double>(histH, "hH");
      }

      Anp::Hist1d<double> histC(*hR);
      Anp::Hist1d<double> histO(*hR);
      
      histC.Reset();
      histO.Reset();

      histC.Fill(*hC);
      histO.Fill(*hO);

      hC = Anp::CreateTH1<double>(histC, "hC");
      hO = Anp::CreateTH1<double>(histO, "hO");
   }

   double miny = -1.0, maxy = -1.0;

   const int rbin = hR -> GetNbinsX();
   
   for(int ibin = 1; ibin <= rbin; ++ibin)
   {      
      const double valueC = hC -> GetBinContent(ibin);
      const double valueO = hO -> GetBinContent(ibin);

      if(!(valueC > 0) || !(valueO > 0))
      {
         continue;
      }

      if(std::fabs(hC ->GetBinLowEdge(ibin)-hO->GetBinLowEdge(ibin)) > 0.00001*hC->GetBinWidth(ibin))
      {
         cerr << "Hist::GetRatio - invalid lower bin edge for hC and hO" << endl;
         continue;
      }

      double errorC = hC -> GetBinError(ibin);

      if(hH && hL)
      {
         if(std::fabs(hC ->GetBinLowEdge(ibin)-hH->GetBinLowEdge(ibin)) > 0.00001*hC->GetBinWidth(ibin) ||
            std::fabs(hC ->GetBinLowEdge(ibin)-hL->GetBinLowEdge(ibin)) > 0.00001*hC->GetBinWidth(ibin))
         {
            cerr << "Hist::GetRatio - invalid lower bin edge for hC and hL or hH" << endl;
            continue;
         }

         const double errorH = hH -> GetBinContent(ibin) - valueC;
         const double errorL = valueC - hL -> GetBinContent(ibin);

         //cout << "errorH = " << errorH << endl;
         //cout << "errorL = " << errorL << endl;

         assert(!(errorH < 0.0) && !(errorL < 0.0) && "Negative errors");

         errorC = std::sqrt(valueC);
         
         if(errorH > errorL)
         {
            errorC += errorH;
         }
         else
         {
            errorC += errorL;
         }
      }

      const double errorO = hO -> GetBinError(ibin);

      const double ratio = valueO/valueC;
      const double error = ratio*(errorO/valueO + errorC/valueC);

      hR -> SetBinContent(ibin, ratio);
      hR -> SetBinError(ibin, error);
   
      if(maxy < 0.0)
      {
         miny = ratio - error;
         maxy = ratio + error;
      }
      else
      {
         miny = std::min<double>(miny, ratio - error);
         maxy = std::max<double>(maxy, ratio + error);
      }
   }

   if(maxy > 0.0)
   {
      const double span = std::fabs(maxy) + std::fabs(miny);
      hR -> GetYaxis() -> SetRangeUser(miny - span*factor, maxy + span*factor);
   }

   if(min_ < max_)
   {
      hR -> GetYaxis() -> SetRangeUser(min_, max_);
   }

   if(option.find("l") != string::npos)
   {
      TLine *line = 0;
      
      if(min_ < max_)
      {
         line = new TLine(min_, 1.0, max_, 1.0);
      }
      else
      {
         line = new TLine(hR -> GetXaxis() -> GetXmin(), 1.0, hR -> GetXaxis() -> GetXmax(), 1.0);
      }

      line -> SetLineWidth(1);
      hR -> GetListOfFunctions() -> Add(line);
   }

   if(option.find("merge") != string::npos)
   {
      if(hO) delete hO;
      if(hC) delete hC;
      if(hL) delete hL;
      if(hH) delete hH;
   }

   return hR;
}

//---------------------------------------------------------------------------------------
std::pair<TH1 *, double> Plot::Hist::Read(const string &hpath)
{
   pair<TH1 *, double> result(0, 0.0);

   if(fFilePath.size() < 6)
   {
      return result;
   }
   
   TFile *file = new TFile(fFilePath.c_str(), "READ");
   if(!file -> IsOpen())
   {
      cerr << "Hist::Read - Failed to open ROOT file " << fFilePath << endl;
      return result;
   }

   TH1 *h = dynamic_cast<TH1 *>(file -> Get(hpath.c_str()));
   if(!h)
   {
      cerr << "Plot::Read - failed to find histogram " << hpath << endl;
      return result;
   }

   h = dynamic_cast<TH1 *>(h -> Clone());
   h -> SetDirectory(0);

   result.first = h;

   DataBlock *block = dynamic_cast<DataBlock *>(file -> Get("block"));   
   if(block)
   {
      result.second = block -> GetProtonsOut();
   }
   else
   {
      cout << "Plot::Get() - failed to find DataBlock in " << file -> GetName() << endl;
   }      

   
   file -> Close();   
   
   return result;
}

//-----------------------------------------------------------------------------------------
TH1* Plot::Hist::CreateRatio(TH1 *hC, TH1 *hO, const string &option, double nmin) const
{
   if(!hC || !hO)
   {
      return 0;
   }

   const int nbin = hC -> GetNbinsX();

   if(hC -> GetNbinsX() != nbin)
   {
      cerr << "Plot::CreateRatio - mismatched number of bins" << endl;
      return 0;
   }

   TH1 *hR = 0;

   if(option.find("merge") != string::npos)
   {
      Anp::Hist1d<double> hist(*hO);

      hist.Merge(nmin);
      hist.Reset();
      hist.Fill(*hC);
      hist.Merge(nmin);

      hR = Anp::CreateTH1<double>(hist, "hR");
   }
   else
   {
      hR = dynamic_cast<TH1 *>(hO -> Clone("hR"));

      if(!hR)
      {
         cerr << "Plot::CreateRatio - failed to clone histogram" << endl;
         return 0;
      }     
   }

   if(hR)
   {
      hR -> Reset();
      hR -> SetDirectory(0);
      hR -> SetStats(false);
      hR -> GetXaxis() -> SetTitle(hO -> GetXaxis() -> GetTitle());
      hR -> GetXaxis() -> CenterTitle();
   }

   return hR;
}

//-----------------------------------------------------------------------------------------
double* Plot::Hist::Create(const int size)
{
   if(size < 1)
   {
      cerr << "Error::Create - can not create zero size array" << endl;
      return 0;
   }

   double *a = new double[size];

   for(int i = 0; i < size; ++i)
   {
      a[i] = 0.0;
   }

   fArray.push_back(a);
   
   return a;
}

//-----------------------------------------------------------------------------------------
void Plot::Hist::ClearError()
{
   if(fAGraph)
   {
      delete fAGraph;
      fAGraph = 0;
   }

   for(vector<double *>::iterator it = fArray.begin(); it != fArray.end(); ++it)
   {
      double *a = *it;
      if(a)
      {
         delete [] a;
      }
   }  

   for(map<int, TBox *>::iterator it = fBox.begin(); it != fBox.end(); ++it)
   {
      TBox *box = it -> second;
      if(box) delete box;
   }

   fArray.clear();
   fError.clear();
   fBox.clear();
}

//-----------------------------------------------------------------------------------------
Plot::Range::Range(const double factor, const double min, const double max)
   :fFactor(factor),
    fMin(min),
    fMax(max)
{
   if(!(fFactor > 1.0))
   {
      fFactor = 1.0;
   }
}

//-----------------------------------------------------------------------------------------
Plot::Range::~Range()
{
}

//-----------------------------------------------------------------------------------------
void Plot::Range::Set(const double factor, const double min, const double max)
{
   fFactor = factor;
   if(min < max)
   {
      fMin = min;
      fMax = max;
   }
}

//-----------------------------------------------------------------------------------------
void Plot::Range::Add(TH1 * h)
{
   if(h)
   {
      fHist.push_back(h);
   }
}

//-----------------------------------------------------------------------------------------
void Plot::Range::Resize(string option, const double min_, const double max_)
{
   if(option.find("Y") != string::npos)
   {
      double min = 0.0, max = 0.0;
      for(HistVec::iterator hit = fHist.begin(); hit != fHist.end(); ++hit)
      {
         TH1 *h = *hit;
         const std::pair<double, double> pmm = GetMinMax(h, option);
         
         if(hit == fHist.begin())
         {
            min = pmm.first;
            max = pmm.second;
         }
         else
         {
            min = std::min(min, pmm.first);
            max = std::max(max, pmm.second);
         }
      }

      if(option.find("min") != string::npos)
      {
         min = min_;
      }
      if(option.find("max") != string::npos)
      {
         max = max_;
      }

      if(min < 0.0)
      {
         min *= fFactor;
      }
      else
      {
         min /= fFactor;
      }
      
      if(max < 0.0)
      {
         max /= fFactor;
      }
      else
      {
         max *= fFactor;
      }
    
      if(option.find("V") != string::npos)
      {
         cout << "Range::Resize - setting y axis range from " << min << " to " << max << endl;
      }

      for(HistVec::iterator hit = fHist.begin(); hit != fHist.end(); ++hit)
      {
         TH1 *h = *hit;
         
         h -> GetYaxis() -> SetRangeUser(min, max);
      }
   }
   else if(option.find("X") != string::npos && min_ < max_)
   {
      if(option.find("V") != string::npos)
      {
         cout << "Range::Resize - setting x axis range from " << min_ << " to " << max_ << endl;
      }

      for(HistVec::iterator hit = fHist.begin(); hit != fHist.end(); ++hit)
      {
         TH1 *h = *hit;

         h -> GetXaxis() -> SetRangeUser(min_, max_);
      }
   }
}


//-----------------------------------------------------------------------------------------
const std::pair<double, double> Plot::Range::GetMinMax(TH1 *h, const string &option) const
{
   if(!h || h -> GetNbinsX() < 1)
   {
      return pair<double, double> (0.0, 0.0);
   }
   
   bool init = false, use_error = false;
   double min = 0.0, max = 0.0;

   if(option.find("E") != string::npos)
   {
      use_error = true;
   }

   for(int ibin = 1; ibin <= h -> GetNbinsX(); ++ibin)
   {
      double value = h -> GetBinContent(ibin);
      double error = h -> GetBinError(ibin);
      
      if(!(value > 0.0))
      {
         continue;
      }
      
      if(use_error)
      {
         value += error;
      }

      if(!init)
      {
         min = value;
         max = value;
         init = true;
         continue;
      }
      
      
      min = std::min(min, value);
      max = std::max(max, value);            
   }

    return pair<double, double> (min, max);
}


//-----------------------------------------------------------------------------------------
Plot::Band::Band()
   :bin(0),
    center(0.0),
    low(0.0),
    high(0.0)
{
}

//-----------------------------------------------------------------------------------------
Plot::Band::Band(int bin_, double center_, double low_, double high_)
   :bin(bin_),
    center(center_),
    low(low_),
    high(high_)
{
}

//-----------------------------------------------------------------------------------------
Plot::Band::~Band()
{
}

//-----------------------------------------------------------------------------------------
int Plot::Band::Bin() const
{
   return bin;
}

//-----------------------------------------------------------------------------------------
double Plot::Band::Center() const
{
   return center;
}

//-----------------------------------------------------------------------------------------
double Plot::Band::Low() const
{
   return low;
}

//-----------------------------------------------------------------------------------------
double Plot::Band::High() const
{
   return high;
}

//-----------------------------------------------------------------------------------------
bool Plot::operator==(const Band &lhs, const Band &rhs)
{
   return (lhs.Bin() == rhs.Bin());
}

//-----------------------------------------------------------------------------------------
void Plot::TestHist(TH1 *hO, TH1 *hC, string option)
{
   const double nmin = 1000;
      
   vector<Anp::Bin<double> > bvec;
   
   map<int, double> bmap;
   
   Anp::Hist1d<double> histO(*hO);
   Anp::Hist1d<double> histC(*hC);
   
   cout << "OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO" << endl;
   histO.Print();
   const unsigned int nmergeO = histO.Merge(1000.0);
   histO.Print();
   
   cout << "Merged " << nmergeO << " bins" << endl;
   
   cout << "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC" << endl;
   histC.Print();
   const unsigned int nmergeC = histC.Merge(1000.0);
   histC.Print();
   
   cout << "Merged " << nmergeC << " bins" << endl;
   
   cout << "NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN" << endl;
   histO.Reset();
   histO.Fill(*hC);
   histO.Print();

   cout << "Attempting draw test..." << endl;

   TH1 *h1 = Anp::CreateTH1<double>(histO, "h1");
   TH1 *h2 = Anp::CreateTH1<double>(histC, "h2");
   
   hC -> SetLineColor(1);
   hC -> SetLineWidth(2);
   hC -> SetFillColor(0);
   h1 -> SetLineColor(2);
   h2 -> SetLineColor(4);

   if(!h1)
   {
      cerr << "Failed to create h1" << endl;
      return;
   }
   if(!h2)
   {
      cerr << "Failed to create h2" << endl;
      return;
   }

   TCanvas *canvas = new TCanvas("ctest", "ctest", 0, 0, 650, 400);
   canvas -> SetGrid();
   canvas -> Draw();
   
   hC -> Draw("HIST");
   h1 -> Draw("sames");   
   h2 -> Draw("sames");

   cout << "Finished draw test" << endl;

   return;

   const int nbin = hC -> GetNbinsX();

   for(int ibin = 1; ibin <= nbin; ++ibin)
   {
      const double value = hO -> GetBinContent(ibin);
      const double edge  = hO -> GetBinLowEdge(ibin);
      
      Anp::Bin<double> bin(ibin, edge);
      bin.add(value);
      
      bvec.push_back(bin);
   }
   
   unsigned int osize = bvec.size();
   double osum = 0.0;
   
   if(option.find("V") != string::npos)
   {
      for(unsigned int ibin = 0; ibin < bvec.size(); ++ibin)
      {
         const Anp::Bin<double> &bin = bvec[ibin];
         cout << "ibin = " << ibin << ", edge = " << bin.edge() << ", sum = " << bin.sum() << endl; 
         osum += bin.sum();
      }
   }

   std::sort(bvec.begin(), bvec.end());
   
   vector<Anp::Bin<double> >::iterator bit = bvec.begin();
   while(bit != bvec.end())
   {
      if(bit -> sum() > nmin || bit == bvec.begin() || bit + 1 == bvec.end())
      {
         ++bit;
         continue;
      }
      
      vector<Anp::Bin<double> >::iterator lit = bit - 1;
      vector<Anp::Bin<double> >::iterator hit = bit + 1;
      
      if(lit -> sum() < hit -> sum())
      {
         Anp::Bin<double> bin(0, lit -> edge());
         
         bin.add(bit -> sum());
         bin.add(lit -> sum());
         
         bvec.erase(lit, hit + 1);
         
         bvec.push_back(bin);
      }
      else
      {
         Anp::Bin<double> bin(0, bit -> edge());
         
         bin.add(bit -> sum());
         bin.add(hit -> sum());
         
         bvec.erase(bit, hit + 1);
         
         bvec.push_back(bin);       
      }
      
      std::sort(bvec.begin(), bvec.end());
      
         bit = bvec.begin();
   }
   
   if(option.find("V") != string::npos)
   {
      double nsum = 0.0;
      for(unsigned int ibin = 0; ibin < bvec.size(); ++ibin)
      {
         const Anp::Bin<double> &bin = bvec[ibin];
         cout << "ibin = " << ibin << ", edge = " << bin.edge() << ", sum = " << bin.sum() << endl; 
         nsum += bin.sum();
      }
      
      cout << "old bvec size = " << osize << endl;
      cout << "new bvec size = " << bvec.size() << endl;
      cout << "osum - nsum = " << osum << " - " << nsum << " = " << osum - nsum << endl;
   }
}

//-----------------------------------------------------------------------------------------
void Plot::TestError(TH1 *hI, string option)
{
   cout << "Plot::TestError - option = " << option << endl;

   TH1 *hO = dynamic_cast<TH1 *>(hI -> Clone("hO"));
   TH1 *hN = dynamic_cast<TH1 *>(hI -> Clone("hN"));

   hO -> SetDirectory(0);
   hN -> SetDirectory(0);

   const int nbin = hN -> GetNbinsX();

   for(int ibin = 1; ibin <= nbin; ++ibin)
   {
      const double value = hN -> GetBinContent(ibin);
      
      if(option.find("half") != std::string::npos)
      {
         hN -> SetBinError(ibin, 0.5*std::sqrt(value));
      }
      else
      {
         hN -> SetBinError(ibin, std::sqrt(value));
      }
   }

   hO -> SetLineColor(2);
   hN -> SetLineColor(4);

   TCanvas *canvas = new TCanvas("ctest", "ctest", 0, 0, 720, 450);
   canvas -> SetGrid();
   canvas -> Draw();
   
   hO -> Draw("E1");
   hN -> Draw("E1 sames");  
}

---