BCH2D Class Reference

#include <BCH2D.h>

List of all members.


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.


Public Member Functions

Constructors and destructors


 BCH2D ()
 ~BCH2D ()
Member functions (miscellaneous methods)


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


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


void SetGlobalMode (double mode[2])
void SetHistogram (TH2D *hist)

Private Attributes

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

Constructor & Destructor Documentation

BCH2D::BCH2D (  ) 

The default constructor.

Definition at line 28 of file BCH2D.cxx.

00029 {
00030    fHistogram = 0;
00031    fIntegratedHistogram = 0;
00032 
00033    fModeFlag = 0;
00034 }

BCH2D::~BCH2D (  ) 

The default destructor.

Definition at line 38 of file BCH2D.cxx.

00039 {
00040    if (fHistogram) delete fHistogram;
00041    if (fIntegratedHistogram) delete fIntegratedHistogram;
00042 }


Member Function Documentation

void BCH2D::CalculateIntegratedHistogram (  ) 

Definition at line 335 of file BCH2D.cxx.

00336 {
00337    int nz = 100;
00338 
00339    double zmax = fHistogram -> GetMaximum();
00340    double dz   = zmax / double(nz);
00341 
00342    double nx = fHistogram -> GetNbinsX();
00343    double ny = fHistogram -> GetNbinsY();
00344 
00345    // create histogram
00346    if (fIntegratedHistogram)
00347       delete fIntegratedHistogram;
00348 
00349    fIntegratedHistogram = new TH1D(
00350          TString::Format("%s_int_prob_%d",fHistogram->GetName(),BCLog::GetHIndex()), "", nz, 0.0, zmax);
00351    fIntegratedHistogram -> SetXTitle("z");
00352    fIntegratedHistogram -> SetYTitle("Integrated probability");
00353    fIntegratedHistogram -> SetStats(kFALSE);
00354 
00355    // loop over histogram
00356    for (int ix = 1; ix <= nx; ix++)
00357       for (int iy = 1; iy <= ny; iy++)
00358       {
00359          int binmin = BCMath::Nint(fHistogram -> GetBinContent(ix, iy) / dz);
00360          for (int i = binmin; i <= nz; i++)
00361             fIntegratedHistogram -> SetBinContent(i,
00362                   fIntegratedHistogram -> GetBinContent(i) +
00363                   fHistogram -> GetBinContent(ix, iy));
00364       }
00365 }

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 82 of file BCH2D.cxx.

00083 {
00084    // draw histogram
00085    fHistogram -> SetLineColor(kBlack);
00086 // fHistogram -> SetLineWidth(4);
00087 
00088    fHistogram->Scale(1./fHistogram->Integral("width"));
00089 
00090    double modex,modey;
00091 
00092    // set mode to display
00093    if(fModeFlag)
00094    {
00095       modex = fMode[0];
00096       modey = fMode[1];
00097    }
00098    else
00099    {
00100       int maximumbin = fHistogram -> GetMaximumBin();
00101 
00102       int binx = maximumbin % (fHistogram -> GetNbinsX() + 2);
00103       int biny = maximumbin / (fHistogram -> GetNbinsX() + 2);
00104 
00105       modex = fHistogram -> GetXaxis() -> GetBinCenter(binx);
00106       modey = fHistogram -> GetYaxis() -> GetBinCenter(biny);
00107    }
00108 
00109    // normalize histogram
00110    fHistogram->Scale(1./fHistogram->Integral("width"));
00111 
00112    // draw according to selected option
00113    if (options == 0)
00114       fHistogram -> Draw("CONT0");
00115    else if (options == 1)
00116    {
00117       fHistogram -> Draw("CONT3");
00118 
00119       // set contours
00120       this -> CalculateIntegratedHistogram();
00121 
00122       double levels[4];
00123       levels[0] = 0.;
00124       levels[1] = this -> GetLevel(1.0 - 0.6827);
00125       levels[2] = this -> GetLevel(1.0 - 0.9545);
00126       levels[3] = this -> GetLevel(1.0 - 0.9973);
00127 
00128       fHistogram -> SetContour(4, levels);
00129 
00130       // best fit value
00131       TMarker* marker = new TMarker(modex, modey, 24);
00132       marker -> Draw();
00133 
00134       TLegend* legend = new TLegend(0.65, 0.80, 0.95, 0.95);
00135       legend -> SetBorderSize(0);
00136       legend -> SetFillColor(kWhite);
00137       legend -> AddEntry(fHistogram, "68% prob. region", "L");
00138       legend -> AddEntry(marker, "Best fit", "P");
00139       legend -> Draw();
00140    }
00141    else if (options == 2)
00142    {
00143       fHistogram -> Draw("CONT3");
00144 
00145       // set contours
00146       this -> CalculateIntegratedHistogram();
00147 
00148       double levels[2];
00149       double level32 = this -> GetLevel(0.32);
00150       levels[0] = 0.;
00151       levels[1] = level32;
00152 
00153       fHistogram -> SetContour(2, levels);
00154 
00155       // best fit value
00156       TMarker* marker = new TMarker(modex, modey, 24);
00157       marker -> Draw();
00158 
00159       TLegend* legend = new TLegend(0.65, 0.80, 0.95, 0.95);
00160       legend -> SetBorderSize(0);
00161       legend -> SetFillColor(kWhite);
00162       legend -> AddEntry(fHistogram, "68% prob. region", "L");
00163       legend -> AddEntry(marker, "Best fit", "P");
00164       legend -> Draw();
00165 
00166    }
00167    else if (options == 3)
00168    {
00169       fHistogram -> Draw("CONT3");
00170 
00171       // set contours
00172       this -> CalculateIntegratedHistogram();
00173 
00174       double levels[2];
00175       double level10 = this -> GetLevel(0.10);
00176       levels[0] = 0.;
00177       levels[1] = level10;
00178 
00179       fHistogram -> SetContour(2, levels);
00180 
00181       TLegend* legend = new TLegend(0.65, 0.80, 0.95, 0.95);
00182       legend -> SetBorderSize(0);
00183       legend -> SetFillColor(kWhite);
00184       legend -> AddEntry(fHistogram, "90% prob. region", "L");
00185       legend -> Draw();
00186    }
00187    else if (options == 4)
00188    {
00189       fHistogram -> Draw("CONT3");
00190 
00191       // set contours
00192       this -> CalculateIntegratedHistogram();
00193 
00194       double levels[2];
00195       double level5 = this -> GetLevel(0.05);
00196       levels[0] = 0.;
00197       levels[1] = level5;
00198 
00199       fHistogram -> SetContour(2, levels);
00200 
00201       TLegend* legend = new TLegend(0.65, 0.80, 0.95, 0.95);
00202       legend -> SetBorderSize(0);
00203       legend -> SetFillColor(kWhite);
00204       legend -> AddEntry(fHistogram, "95% prob. region", "L");
00205       legend -> Draw();
00206    }
00207    else if (options == 5)
00208       fHistogram -> Draw("COL");
00209    else if (options == 52 || options == 521)
00210    {
00211       // create new empty histogram
00212       int nx = fHistogram -> GetNbinsX();
00213       int ny = fHistogram -> GetNbinsY();
00214       double xmin = fHistogram->GetXaxis()->GetXmin();
00215       double xmax = fHistogram->GetXaxis()->GetXmax();
00216       TH2D * h = new TH2D(
00217             TString::Format("htemp52_%d",BCLog::GetHIndex()),fHistogram->GetTitle(),
00218             nx,xmin,xmax,
00219             ny,fHistogram->GetYaxis()->GetXmin(),fHistogram->GetYaxis()->GetXmax());
00220       h->SetXTitle(fHistogram->GetXaxis()->GetTitle());
00221       h->SetYTitle(fHistogram->GetYaxis()->GetTitle());
00222 
00223       // copy contents of the main histogram
00224 //    double min = fHistogram->GetMinimum(0.);
00225       for(int i=1;i<=nx;i++)
00226          for(int j=1;j<=ny;j++)
00227          {
00228             double val = fHistogram -> GetBinContent(i,j);
00229             // if requested, change contents of bins to log scale
00230             if(options == 521)
00231             {
00232 //             if(val == 0.)
00233 //                val = log(min)-1.;
00234 //             else
00235                   val = log10(val);
00236             }
00237             h -> SetBinContent(i,j,val);
00238          }
00239 
00240       // draw
00241       h->SetStats(0);
00242       h->Draw("col");
00243 
00244       // draw contour
00245       fHistogram -> Draw("cont3 same");
00246       fHistogram -> SetLineWidth(2);
00247 
00248       // set contours
00249       this -> CalculateIntegratedHistogram();
00250 
00251       double levels[] = { this -> GetLevel(0.32) };
00252       fHistogram -> SetContour(1, levels);
00253 
00254       // best fit value
00255       if(drawmode)
00256       {
00257          TMarker * marker0 = new TMarker(modex, modey, 8);
00258          marker0 -> SetMarkerColor(0);
00259          marker0 -> SetMarkerSize(.7);
00260          marker0 -> Draw();
00261          TMarker * marker1 = new TMarker(modex, modey, 4);
00262          marker1 -> SetMarkerColor(1);
00263          marker1 -> SetMarkerSize(.7);
00264          marker1 -> Draw();
00265 //       TMarker * marker = new TMarker(modex, modey, 5);
00266 //       marker -> SetMarkerColor(0);
00267 //       marker -> Draw();
00268       }
00269    }
00270    else if (options == 53 || options == 531)
00271    {
00272       // create new empty histogram
00273 //    int nx = fHistogram -> GetNbinsX();
00274 //    int ny = fHistogram -> GetNbinsY();
00275 //    double xmin = fHistogram->GetXaxis()->GetXmin();
00276 //    double xmax = fHistogram->GetXaxis()->GetXmax();
00277 //    TH2D * h = new TH2D(
00278 //          TString::Format("htemp53_%d",BCLog::GetHIndex()),fHistogram->GetTitle(),
00279 //          nx,xmin,xmax,
00280 //          ny,fHistogram->GetYaxis()->GetXmin(),fHistogram->GetYaxis()->GetXmax());
00281 //    h->SetXTitle(fHistogram->GetXaxis()->GetTitle());
00282 //    h->SetYTitle(fHistogram->GetYaxis()->GetTitle());
00283 
00284       // copy contents of the main histogram
00285 //    double min = fHistogram->GetMinimum(0.);
00286 /*    for(int i=1;i<=nx;i++)
00287          for(int j=1;j<=ny;j++)
00288          {
00289             double val = fHistogram -> GetBinContent(i,j);
00290             // if requested, change contents of bins to log scale
00291             if(options == 531)
00292             {
00293 //             if(val == 0.)
00294 //                val = log(min)-1.;
00295 //             else
00296                   val = log10(val);
00297             }
00298             h -> SetBinContent(i,j,val);
00299          }
00300 
00301       // draw
00302       h->SetStats(0);
00303       h->Draw("colz");
00304 */
00305       gPad -> SetLogz();
00306       fHistogram -> Draw("colz");
00307 
00308       // draw contour
00309 //    fHistogram -> Draw("cont3 same");
00310 //    fHistogram -> SetLineWidth(2);
00311 
00312       // set contours
00313 //    this -> CalculateIntegratedHistogram();
00314 
00315 //    double levels[] = { this -> GetLevel(0.32) };
00316 //    fHistogram -> SetContour(1, levels);
00317 
00318       // best fit value
00319       TMarker * marker0 = new TMarker(modex, modey, 8);
00320       marker0 -> SetMarkerColor(0);
00321       marker0 -> SetMarkerSize(.7);
00322       marker0 -> Draw();
00323       TMarker * marker1 = new TMarker(modex, modey, 4);
00324       marker1 -> SetMarkerColor(1);
00325       marker1 -> SetMarkerSize(.7);
00326       marker1 -> Draw();
00327 //    TMarker * marker = new TMarker(modex, modey, 5);
00328 //    marker -> SetMarkerColor(0);
00329 //    marker -> Draw();
00330    }
00331 }

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

Definition at line 538 of file BCH2D.cxx.

00539 {
00540    // define new graph
00541    int nx = h -> GetNbinsX() - 1;
00542 
00543    TGraph * g = new TGraph(2*nx);
00544 // g -> SetFillStyle(1001);
00545 // g -> SetFillColor(kYellow);
00546 
00547    // get error bands
00548    std::vector <double> ymin = GetLevelBoundary(h,l1);
00549    std::vector <double> ymax = GetLevelBoundary(h,l2);
00550 
00551    for (int i = 0; i < nx; i++)
00552    {
00553       g -> SetPoint(i, h -> GetXaxis() -> GetBinCenter(i+1), ymin[i]);
00554       g -> SetPoint(nx+i, h -> GetXaxis() -> GetBinCenter(nx-i), ymax[nx-i-1]);
00555    }
00556 
00557    return g;
00558 }

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

Definition at line 382 of file BCH2D.cxx.

00383 {
00384    n=0;
00385 
00386    int nbands=0;
00387    TH2D * hcopy = (TH2D*)h->Clone(TString::Format("%s_copy_%d",h->GetName(),BCLog::GetHIndex()));
00388 
00389    std::vector <int> nint=GetNIntervalsY(hcopy,nbands);
00390 
00391    if(nbands>2)
00392    {
00393       BCLog::Out(BCLog::warning,BCLog::warning,
00394          Form("BCH2D::GetBandGraphs : Detected %d bands. Maximum allowed is 2 (sorry).",nbands));
00395       return 0;
00396    }
00397    else if(nbands==0)
00398    {
00399       BCLog::Out(BCLog::warning,BCLog::warning,"BCH2D::GetBandGraphs : No bands detected.");
00400       return 0;
00401    }
00402 
00403    TGraph ** gxx = new TGraph*[nbands];
00404 
00405    TH2D * h0 = (TH2D*)h->Clone();
00406 
00407    if (nbands>0)
00408       gxx[0] = GetLowestBandGraph(h0,nint);
00409 
00410    if (nbands==2)
00411    {
00412       gxx[1] = GetLowestBandGraph(h0);
00413       n=2;
00414    }
00415    else
00416       n=1;
00417 
00418    return gxx;
00419 }

TH2D* BCH2D::GetHistogram (  )  [inline]

Returns:
The 2D histogram.

Definition at line 55 of file BCH2D.h.

00056          { return fHistogram; };

double BCH2D::GetLevel ( double  p  ) 

Definition at line 369 of file BCH2D.cxx.

00370 {
00371    double quantiles[1];
00372    double probsum[1];
00373    probsum[0] = p;
00374 
00375    fIntegratedHistogram -> GetQuantiles( 1, quantiles, probsum);
00376 
00377    return quantiles[0];
00378 }

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

Definition at line 511 of file BCH2D.cxx.

00512 {
00513    std::vector <double> b;
00514 
00515    int nx = h -> GetNbinsX();
00516 
00517    b.assign(nx - 1, 0.0);
00518 
00519    // loop over x and y bins.
00520    for (int ix = 1; ix < nx; ix++)
00521    {
00522       TH1D * h1 = h -> ProjectionY(TString::Format("temphist_%d",BCLog::GetHIndex()), ix, ix);
00523 
00524       int nprobSum = 1;
00525       double q[1];
00526       double probSum[] = { level };
00527 
00528       h1 -> GetQuantiles(nprobSum, q, probSum);
00529 
00530       b[ix-1] = q[0];
00531    }
00532 
00533    return b;
00534 }

TGraph * BCH2D::GetLowestBandGraph ( TH2D *  h  ) 

Definition at line 459 of file BCH2D.cxx.

00460 {
00461    int n;
00462    return GetLowestBandGraph(h,GetNIntervalsY(h,n));
00463 }

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

Definition at line 467 of file BCH2D.cxx.

00468 {
00469    int nx = h -> GetNbinsX() - 1;
00470    int ny = h -> GetNbinsY();
00471 
00472    TGraph * g = new TGraph(2*nx);
00473 
00474    for (int ix=1; ix<=nx; ix++)
00475    {
00476       // get x for the bin
00477       double x;
00478       if(ix==1)
00479          x = h->GetXaxis()->GetBinLowEdge(1);
00480       else if(ix==nx)
00481          x = h->GetXaxis()->GetBinLowEdge(nx+1);
00482       else
00483          x = h->GetXaxis()->GetBinCenter(ix);
00484 
00485       for(int iy=1; iy<=ny; iy++)
00486          if(h->GetBinContent(ix,iy)>0.)
00487          {
00488             // get low edge of the first not empty bin in y
00489             g->SetPoint(ix-1, x, h->GetYaxis()->GetBinLowEdge(iy));
00490 
00491             // delete content of all subsequent not empty bins
00492             if(nint[ix-1]==2)
00493                h->SetBinContent(ix,iy,0.);
00494 
00495             while(h->GetBinContent(ix,++iy)>0.)
00496                if(nint[ix-1]==2)
00497                   h->SetBinContent(ix,iy,0.);
00498 
00499             // get low edge of the first empty bin in y
00500             g->SetPoint(2*nx-ix, x, h->GetYaxis()->GetBinLowEdge(iy));
00501 
00502             break;
00503          }
00504    }
00505 
00506    return g;
00507 }

void BCH2D::GetMode ( double &  mode  ) 

Returns:
The mean of the distribution.

The mode of the distribution

Definition at line 46 of file BCH2D.cxx.

00047 {
00048    // what is this?
00049 
00050 // double mode[2];
00051 // int binx, biny, binz;
00052 
00053 // fHistogram -> GetMaximumBin(binx, biny, binz);
00054 // mode = fHistogram -> GetBinCenter();
00055 
00056 // return mode;
00057 }

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 423 of file BCH2D.cxx.

00424 {
00425    std::vector <int> nint;
00426 
00427    int nx = h -> GetNbinsX();
00428    int ny = h -> GetNbinsY();
00429 
00430    nfoundmax=0;
00431 
00432    // loop over histogram bins in x
00433    for (int ix=1; ix<=nx; ix++)
00434    {
00435       int nfound=0;
00436 
00437       // loop over histogram bins in y
00438       // count nonzero intervals in y
00439       for (int iy=1; iy<=ny; iy++)
00440          if(h->GetBinContent(ix,iy)>0.)
00441          {
00442             while(h->GetBinContent(ix,++iy)>0.)
00443                ;
00444             nfound++;
00445          }
00446 
00447       // store maximum number of nonzero intervals for the histogram
00448       if(nfound>nfoundmax)
00449          nfoundmax=nfound;
00450 
00451       nint.push_back(nfound);
00452    }
00453 
00454    return nint;
00455 }

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 61 of file BCH2D.cxx.

00062 {
00063    // create temporary canvas
00064    TCanvas * c;
00065    if(ww >0 && wh > 0)
00066       c = new TCanvas("c","c",ww,wh);
00067    else
00068       c = new TCanvas("c","c",700,700);
00069    c -> cd();
00070 
00071    // draw histogram
00072    this->Draw(options);
00073 
00074    gPad->RedrawAxis();
00075 
00076    // print to file
00077    c -> Print(filename);
00078 }

void BCH2D::SetGlobalMode ( double  mode[2]  )  [inline]

Set global mode.

Parameters:
The global mode.

Definition at line 78 of file BCH2D.h.

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

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

Set the 2D histogram.

Definition at line 72 of file BCH2D.h.

00073          { fHistogram = hist; };


Member Data Documentation

TH2D* BCH2D::fHistogram [private]

The 2D histogram

Definition at line 136 of file BCH2D.h.

TH1D* BCH2D::fIntegratedHistogram [private]

The integrated 2D histogram

Definition at line 140 of file BCH2D.h.

double BCH2D::fMode[2] [private]

Global mode

Definition at line 144 of file BCH2D.h.

int BCH2D::fModeFlag [private]

"Is there a global mode?" flag

Definition at line 148 of file BCH2D.h.


Generated on Fri Jan 16 10:24:10 2009 for BAT - Bayesian Analysis Toolkit by  doxygen 1.5.6