BCH2D.cxx

Go to the documentation of this file.

/*
 * Copyright (C) 2008, Daniel Kollar and Kevin Kroeninger.
 * All rights reserved.
 *
 * For the licensing terms see doc/COPYING.
 */

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

#include "BAT/BCH2D.h"

#include "BAT/BCMath.h"
#include "BAT/BCLog.h"

#include <TH1D.h>
#include <TH2D.h>
#include <TGraph.h>
#include <TCanvas.h>
#include <TLine.h>
#include <TMarker.h>
#include <TLegend.h>
#include <TString.h>

#include <math.h>

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

BCH2D::BCH2D()
{
   fHistogram = 0;
   fIntegratedHistogram = 0;

   fModeFlag = 0;
}

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

BCH2D::~BCH2D()
{
   if (fHistogram) delete fHistogram;
   if (fIntegratedHistogram) delete fIntegratedHistogram;
}

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

void BCH2D::GetMode(double& mode)
{
   // what is this?

// double mode[2];
// int binx, biny, binz;

// fHistogram -> GetMaximumBin(binx, biny, binz);
// mode = fHistogram -> GetBinCenter();

// return mode;
}

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

void BCH2D::Print(const char * filename, int options, int ww, int wh)
{
   // create temporary canvas
   TCanvas * c;
   if(ww >0 && wh > 0)
      c = new TCanvas("c","c",ww,wh);
   else
      c = new TCanvas("c","c",700,700);
   c -> cd();

   // draw histogram
   this->Draw(options);

   gPad->RedrawAxis();

   // print to file
   c -> Print(filename);
}

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

void BCH2D::Draw(int options, bool drawmode)
{
   // draw histogram
   fHistogram -> SetLineColor(kBlack);
// fHistogram -> SetLineWidth(4);

   fHistogram->Scale(1./fHistogram->Integral("width"));

   double modex,modey;

   // set mode to display
   if(fModeFlag)
   {
      modex = fMode[0];
      modey = fMode[1];
   }
   else
   {
      int maximumbin = fHistogram -> GetMaximumBin();

      int binx = maximumbin % (fHistogram -> GetNbinsX() + 2);
      int biny = maximumbin / (fHistogram -> GetNbinsX() + 2);

      modex = fHistogram -> GetXaxis() -> GetBinCenter(binx);
      modey = fHistogram -> GetYaxis() -> GetBinCenter(biny);
   }

   // normalize histogram
   fHistogram->Scale(1./fHistogram->Integral("width"));

   // draw according to selected option
   if (options == 0)
      fHistogram -> Draw("CONT0");
   else if (options == 1)
   {
      fHistogram -> Draw("CONT3");

      // set contours
      this -> CalculateIntegratedHistogram();

      double levels[4];
      levels[0] = 0.;
      levels[1] = this -> GetLevel(1.0 - 0.6827);
      levels[2] = this -> GetLevel(1.0 - 0.9545);
      levels[3] = this -> GetLevel(1.0 - 0.9973);

      fHistogram -> SetContour(4, levels);

      // best fit value
      TMarker* marker = new TMarker(modex, modey, 24);
      marker -> Draw();

      TLegend* legend = new TLegend(0.65, 0.80, 0.95, 0.95);
      legend -> SetBorderSize(0);
      legend -> SetFillColor(kWhite);
      legend -> AddEntry(fHistogram, "68% prob. region", "L");
      legend -> AddEntry(marker, "Best fit", "P");
      legend -> Draw();
   }
   else if (options == 2)
   {
      fHistogram -> Draw("CONT3");

      // set contours
      this -> CalculateIntegratedHistogram();

      double levels[2];
      double level32 = this -> GetLevel(0.32);
      levels[0] = 0.;
      levels[1] = level32;

      fHistogram -> SetContour(2, levels);

      // best fit value
      TMarker* marker = new TMarker(modex, modey, 24);
      marker -> Draw();

      TLegend* legend = new TLegend(0.65, 0.80, 0.95, 0.95);
      legend -> SetBorderSize(0);
      legend -> SetFillColor(kWhite);
      legend -> AddEntry(fHistogram, "68% prob. region", "L");
      legend -> AddEntry(marker, "Best fit", "P");
      legend -> Draw();

   }
   else if (options == 3)
   {
      fHistogram -> Draw("CONT3");

      // set contours
      this -> CalculateIntegratedHistogram();

      double levels[2];
      double level10 = this -> GetLevel(0.10);
      levels[0] = 0.;
      levels[1] = level10;

      fHistogram -> SetContour(2, levels);

      TLegend* legend = new TLegend(0.65, 0.80, 0.95, 0.95);
      legend -> SetBorderSize(0);
      legend -> SetFillColor(kWhite);
      legend -> AddEntry(fHistogram, "90% prob. region", "L");
      legend -> Draw();
   }
   else if (options == 4)
   {
      fHistogram -> Draw("CONT3");

      // set contours
      this -> CalculateIntegratedHistogram();

      double levels[2];
      double level5 = this -> GetLevel(0.05);
      levels[0] = 0.;
      levels[1] = level5;

      fHistogram -> SetContour(2, levels);

      TLegend* legend = new TLegend(0.65, 0.80, 0.95, 0.95);
      legend -> SetBorderSize(0);
      legend -> SetFillColor(kWhite);
      legend -> AddEntry(fHistogram, "95% prob. region", "L");
      legend -> Draw();
   }
   else if (options == 5)
      fHistogram -> Draw("COL");
   else if (options == 52 || options == 521)
   {
      // create new empty histogram
      int nx = fHistogram -> GetNbinsX();
      int ny = fHistogram -> GetNbinsY();
      double xmin = fHistogram->GetXaxis()->GetXmin();
      double xmax = fHistogram->GetXaxis()->GetXmax();
      TH2D * h = new TH2D(
            TString::Format("htemp52_%d",BCLog::GetHIndex()),fHistogram->GetTitle(),
            nx,xmin,xmax,
            ny,fHistogram->GetYaxis()->GetXmin(),fHistogram->GetYaxis()->GetXmax());
      h->SetXTitle(fHistogram->GetXaxis()->GetTitle());
      h->SetYTitle(fHistogram->GetYaxis()->GetTitle());

      // copy contents of the main histogram
//    double min = fHistogram->GetMinimum(0.);
      for(int i=1;i<=nx;i++)
         for(int j=1;j<=ny;j++)
         {
            double val = fHistogram -> GetBinContent(i,j);
            // if requested, change contents of bins to log scale
            if(options == 521)
            {
//             if(val == 0.)
//                val = log(min)-1.;
//             else
                  val = log10(val);
            }
            h -> SetBinContent(i,j,val);
         }

      // draw
      h->SetStats(0);
      h->Draw("col");

      // draw contour
      fHistogram -> Draw("cont3 same");
      fHistogram -> SetLineWidth(2);

      // set contours
      this -> CalculateIntegratedHistogram();

      double levels[] = { this -> GetLevel(0.32) };
      fHistogram -> SetContour(1, levels);

      // best fit value
      if(drawmode)
      {
         TMarker * marker0 = new TMarker(modex, modey, 8);
         marker0 -> SetMarkerColor(0);
         marker0 -> SetMarkerSize(.7);
         marker0 -> Draw();
         TMarker * marker1 = new TMarker(modex, modey, 4);
         marker1 -> SetMarkerColor(1);
         marker1 -> SetMarkerSize(.7);
         marker1 -> Draw();
//       TMarker * marker = new TMarker(modex, modey, 5);
//       marker -> SetMarkerColor(0);
//       marker -> Draw();
      }
   }
   else if (options == 53 || options == 531)
   {
      // create new empty histogram
//    int nx = fHistogram -> GetNbinsX();
//    int ny = fHistogram -> GetNbinsY();
//    double xmin = fHistogram->GetXaxis()->GetXmin();
//    double xmax = fHistogram->GetXaxis()->GetXmax();
//    TH2D * h = new TH2D(
//          TString::Format("htemp53_%d",BCLog::GetHIndex()),fHistogram->GetTitle(),
//          nx,xmin,xmax,
//          ny,fHistogram->GetYaxis()->GetXmin(),fHistogram->GetYaxis()->GetXmax());
//    h->SetXTitle(fHistogram->GetXaxis()->GetTitle());
//    h->SetYTitle(fHistogram->GetYaxis()->GetTitle());

      // copy contents of the main histogram
//    double min = fHistogram->GetMinimum(0.);
/*    for(int i=1;i<=nx;i++)
         for(int j=1;j<=ny;j++)
         {
            double val = fHistogram -> GetBinContent(i,j);
            // if requested, change contents of bins to log scale
            if(options == 531)
            {
//             if(val == 0.)
//                val = log(min)-1.;
//             else
                  val = log10(val);
            }
            h -> SetBinContent(i,j,val);
         }

      // draw
      h->SetStats(0);
      h->Draw("colz");
*/
      gPad -> SetLogz();
      fHistogram -> Draw("colz");

      // draw contour
//    fHistogram -> Draw("cont3 same");
//    fHistogram -> SetLineWidth(2);

      // set contours
//    this -> CalculateIntegratedHistogram();

//    double levels[] = { this -> GetLevel(0.32) };
//    fHistogram -> SetContour(1, levels);

      // best fit value
      TMarker * marker0 = new TMarker(modex, modey, 8);
      marker0 -> SetMarkerColor(0);
      marker0 -> SetMarkerSize(.7);
      marker0 -> Draw();
      TMarker * marker1 = new TMarker(modex, modey, 4);
      marker1 -> SetMarkerColor(1);
      marker1 -> SetMarkerSize(.7);
      marker1 -> Draw();
//    TMarker * marker = new TMarker(modex, modey, 5);
//    marker -> SetMarkerColor(0);
//    marker -> Draw();
   }
}

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

void BCH2D::CalculateIntegratedHistogram()
{
   int nz = 100;

   double zmax = fHistogram -> GetMaximum();
   double dz   = zmax / double(nz);

   double nx = fHistogram -> GetNbinsX();
   double ny = fHistogram -> GetNbinsY();

   // create histogram
   if (fIntegratedHistogram)
      delete fIntegratedHistogram;

   fIntegratedHistogram = new TH1D(
         TString::Format("%s_int_prob_%d",fHistogram->GetName(),BCLog::GetHIndex()), "", nz, 0.0, zmax);
   fIntegratedHistogram -> SetXTitle("z");
   fIntegratedHistogram -> SetYTitle("Integrated probability");
   fIntegratedHistogram -> SetStats(kFALSE);

   // loop over histogram
   for (int ix = 1; ix <= nx; ix++)
      for (int iy = 1; iy <= ny; iy++)
      {
         int binmin = BCMath::Nint(fHistogram -> GetBinContent(ix, iy) / dz);
         for (int i = binmin; i <= nz; i++)
            fIntegratedHistogram -> SetBinContent(i,
                  fIntegratedHistogram -> GetBinContent(i) +
                  fHistogram -> GetBinContent(ix, iy));
      }
}

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

double BCH2D::GetLevel(double p)
{
   double quantiles[1];
   double probsum[1];
   probsum[0] = p;

   fIntegratedHistogram -> GetQuantiles( 1, quantiles, probsum);

   return quantiles[0];
}

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

TGraph ** BCH2D::GetBandGraphs(TH2D * h, int &n)
{
   n=0;

   int nbands=0;
   TH2D * hcopy = (TH2D*)h->Clone(TString::Format("%s_copy_%d",h->GetName(),BCLog::GetHIndex()));

   std::vector <int> nint=GetNIntervalsY(hcopy,nbands);

   if(nbands>2)
   {
      BCLog::Out(BCLog::warning,BCLog::warning,
         Form("BCH2D::GetBandGraphs : Detected %d bands. Maximum allowed is 2 (sorry).",nbands));
      return 0;
   }
   else if(nbands==0)
   {
      BCLog::Out(BCLog::warning,BCLog::warning,"BCH2D::GetBandGraphs : No bands detected.");
      return 0;
   }

   TGraph ** gxx = new TGraph*[nbands];

   TH2D * h0 = (TH2D*)h->Clone();

   if (nbands>0)
      gxx[0] = GetLowestBandGraph(h0,nint);

   if (nbands==2)
   {
      gxx[1] = GetLowestBandGraph(h0);
      n=2;
   }
   else
      n=1;

   return gxx;
}

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

std::vector <int> BCH2D::GetNIntervalsY(TH2D * h, int &nfoundmax)
{
   std::vector <int> nint;

   int nx = h -> GetNbinsX();
   int ny = h -> GetNbinsY();

   nfoundmax=0;

   // loop over histogram bins in x
   for (int ix=1; ix<=nx; ix++)
   {
      int nfound=0;

      // loop over histogram bins in y
      // count nonzero intervals in y
      for (int iy=1; iy<=ny; iy++)
         if(h->GetBinContent(ix,iy)>0.)
         {
            while(h->GetBinContent(ix,++iy)>0.)
               ;
            nfound++;
         }

      // store maximum number of nonzero intervals for the histogram
      if(nfound>nfoundmax)
         nfoundmax=nfound;

      nint.push_back(nfound);
   }

   return nint;
}

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

TGraph * BCH2D::GetLowestBandGraph(TH2D * h)
{
   int n;
   return GetLowestBandGraph(h,GetNIntervalsY(h,n));
}

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

TGraph * BCH2D::GetLowestBandGraph(TH2D * h, std::vector<int> nint)
{
   int nx = h -> GetNbinsX() - 1;
   int ny = h -> GetNbinsY();

   TGraph * g = new TGraph(2*nx);

   for (int ix=1; ix<=nx; ix++)
   {
      // get x for the bin
      double x;
      if(ix==1)
         x = h->GetXaxis()->GetBinLowEdge(1);
      else if(ix==nx)
         x = h->GetXaxis()->GetBinLowEdge(nx+1);
      else
         x = h->GetXaxis()->GetBinCenter(ix);

      for(int iy=1; iy<=ny; iy++)
         if(h->GetBinContent(ix,iy)>0.)
         {
            // get low edge of the first not empty bin in y
            g->SetPoint(ix-1, x, h->GetYaxis()->GetBinLowEdge(iy));

            // delete content of all subsequent not empty bins
            if(nint[ix-1]==2)
               h->SetBinContent(ix,iy,0.);

            while(h->GetBinContent(ix,++iy)>0.)
               if(nint[ix-1]==2)
                  h->SetBinContent(ix,iy,0.);

            // get low edge of the first empty bin in y
            g->SetPoint(2*nx-ix, x, h->GetYaxis()->GetBinLowEdge(iy));

            break;
         }
   }

   return g;
}

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

std::vector <double> BCH2D::GetLevelBoundary(TH2D * h, double level)
{
   std::vector <double> b;

   int nx = h -> GetNbinsX();

   b.assign(nx - 1, 0.0);

   // loop over x and y bins.
   for (int ix = 1; ix < nx; ix++)
   {
      TH1D * h1 = h -> ProjectionY(TString::Format("temphist_%d",BCLog::GetHIndex()), ix, ix);

      int nprobSum = 1;
      double q[1];
      double probSum[] = { level };

      h1 -> GetQuantiles(nprobSum, q, probSum);

      b[ix-1] = q[0];
   }

   return b;
}

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

TGraph * BCH2D::GetBandGraph(TH2D * h , double l1, double l2)
{
   // define new graph
   int nx = h -> GetNbinsX() - 1;

   TGraph * g = new TGraph(2*nx);
// g -> SetFillStyle(1001);
// g -> SetFillColor(kYellow);

   // get error bands
   std::vector <double> ymin = GetLevelBoundary(h,l1);
   std::vector <double> ymax = GetLevelBoundary(h,l2);

   for (int i = 0; i < nx; i++)
   {
      g -> SetPoint(i, h -> GetXaxis() -> GetBinCenter(i+1), ymin[i]);
      g -> SetPoint(nx+i, h -> GetXaxis() -> GetBinCenter(nx-i), ymax[nx-i-1]);
   }

   return g;
}

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

---