Private Attributes

BCH2D Class Reference

A class for handling 2D distributions. More...

#include <BCH2D.h>

List of all members.

Public Member Functions

Constructors and destructors

 BCH2D ()
 BCH2D (TH2D *h)
 ~BCH2D ()
Member functions (get)

TH2D * GetHistogram ()
void GetMode (double &mode)
Member functions (set)

void SetHistogram (TH2D *hist)
void SetGlobalMode (double mode[2])
Member functions (miscellaneous methods)

void Print (const char *filename, int options=0, int ww=0, int wh=0)
void Draw (int options=0, bool drawmode=true)
void CalculateIntegratedHistogram ()
double GetLevel (double p)
std::vector< int > GetNIntervalsY (TH2D *h, int &nfoundmax)
TGraph * GetLowestBandGraph (TH2D *h, std::vector< int > nint)
TGraph * GetLowestBandGraph (TH2D *h)
std::vector< double > GetLevelBoundary (double level)
std::vector< double > GetLevelBoundary (TH2D *h, double level)
TGraph * GetBandGraph (double level1, double level2)
TGraph * GetBandGraph (TH2D *h, double level1, double level2)
TGraph ** GetBandGraphs (TH2D *h, int &n)

Private Attributes

TH2D * fHistogram
TH1D * fIntegratedHistogram
double fMode [2]
int fModeFlag

Detailed Description

A class for handling 2D distributions.

Author:
Daniel Kollar
Kevin Kröninger
Version:
1.0
Date:
08.2008 This class contains a TH2D histogram and some additional functions. It is used for marginalized distributions.

Definition at line 33 of file BCH2D.h.

Constructor & Destructor Documentation

BCH2D::BCH2D (  ) 

The default constructor.

Definition at line 28 of file BCH2D.cxx.

{
   fHistogram = 0;
   fIntegratedHistogram = 0;

   fModeFlag = 0;
}

BCH2D::BCH2D ( TH2D *  h  ) 

The complete constructor.

Definition at line 38 of file BCH2D.cxx.

{
   fHistogram = h;
   fIntegratedHistogram = 0;

   fModeFlag = 0;
}

BCH2D::~BCH2D (  ) 

The default destructor.

Definition at line 48 of file BCH2D.cxx.

{
   if (fHistogram) delete fHistogram;
   if (fIntegratedHistogram) delete fIntegratedHistogram;
}

Member Function Documentation

void BCH2D::CalculateIntegratedHistogram (  ) 

Definition at line 345 of file BCH2D.cxx.

{
   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));
      }
}

void BCH2D::Draw ( int  options = 0,
bool  drawmode = true 
)

Draw 2-d distribution into the active canvas

Parameters:
options explanation to come
drawmode specify whether a marker should be drawn at the location of the mode

Definition at line 92 of file BCH2D.cxx.

{
   // 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();
   }
}

TGraph * BCH2D::GetBandGraph ( double  level1,
double  level2 
)

Definition at line 555 of file BCH2D.cxx.

{
   return GetBandGraph(fHistogram , l1, l2);
}

TGraph * BCH2D::GetBandGraph ( TH2D *  h,
double  level1,
double  level2 
)

Definition at line 562 of file BCH2D.cxx.

{
   // 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;
}

TGraph ** BCH2D::GetBandGraphs ( TH2D *  h,
int &  n 
)

Definition at line 392 of file BCH2D.cxx.

{
   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;
}

TH2D* BCH2D::GetHistogram (  )  [inline]
Returns:
The 2D histogram.

Definition at line 59 of file BCH2D.h.

         { return fHistogram; };

double BCH2D::GetLevel ( double  p  ) 

Definition at line 379 of file BCH2D.cxx.

{
   double quantiles[1];
   double probsum[1];
   probsum[0] = p;

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

   return quantiles[0];
}

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

Definition at line 528 of file BCH2D.cxx.

{
   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;
}

std::vector< double > BCH2D::GetLevelBoundary ( double  level  ) 

Definition at line 521 of file BCH2D.cxx.

{
   return GetLevelBoundary(fHistogram, level);
}

TGraph * BCH2D::GetLowestBandGraph ( TH2D *  h  ) 

Definition at line 469 of file BCH2D.cxx.

{
   int n;
   return GetLowestBandGraph(h,GetNIntervalsY(h,n));
}

TGraph * BCH2D::GetLowestBandGraph ( TH2D *  h,
std::vector< int >  nint 
)

Definition at line 477 of file BCH2D.cxx.

{
   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;
}

void BCH2D::GetMode ( double &  mode  ) 
Returns:
The mean of the distribution.
The mode of the distribution

Definition at line 56 of file BCH2D.cxx.

{
   // what is this?

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

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

// return mode;
}

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

Returns the number of intervals as a function of x

Parameters:
h The histogram.
nfoundmax The maximum number of intervals.
Returns:
A vector containing the number of intervals for all bins in x.

Definition at line 433 of file BCH2D.cxx.

{
   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;
}

void BCH2D::Print ( const char *  filename,
int  options = 0,
int  ww = 0,
int  wh = 0 
)

Print 2-d histogram to file

Parameters:
filename The filename
ww canvas size in pixels along X
ww canvas size in pixels along Y If ww and wh are set to 0, default ROOT canvas size is used. For explanation of the parameter options see the Draw() method.

Definition at line 71 of file BCH2D.cxx.

{
   // 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::SetGlobalMode ( double  mode[2]  )  [inline]

Set global mode.

Parameters:
The global mode.

Definition at line 82 of file BCH2D.h.

         { fMode[0] = mode[0]; fMode[1] = mode[1]; fModeFlag =1; };

void BCH2D::SetHistogram ( TH2D *  hist  )  [inline]

Set the 2D histogram.

Definition at line 76 of file BCH2D.h.

         { fHistogram = hist; };

Member Data Documentation

TH2D* BCH2D::fHistogram [private]

The 2D histogram

Definition at line 142 of file BCH2D.h.

TH1D* BCH2D::fIntegratedHistogram [private]

The integrated 2D histogram

Definition at line 146 of file BCH2D.h.

double BCH2D::fMode[2] [private]

Global mode

Definition at line 150 of file BCH2D.h.

int BCH2D::fModeFlag [private]

"Is there a global mode?" flag

Definition at line 154 of file BCH2D.h.

The documentation for this class was generated from the following files:
Generated on Mon Aug 30 2010 22:14:55 for Bayesian Analysis Toolkit by  doxygen 1.7.1