BCH1D Class Reference

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

#include <BCH1D.h>

List of all members.

Public Member Functions

Constructors and destructors



 BCH1D (TH1D *hist)
 BCH1D ()
 ~BCH1D ()
Member functions (miscellaneous methods)



void Draw (int options=0, double ovalue=0.)
void DrawShadedLimits (double mode, double min, double max, double limit=0)
void DrawSmallest (double mode, double prob, bool drawmean=true)
double GetSmallestInterval (double &min, double &max, double content=.68)
TH1D * GetSmallestIntervalHistogram (double level)
std::vector< double > GetSmallestIntervals (double content=0.68)
TH1D * GetSubHisto (double min, double max, const char *name)
double IntegralWidth (double min, double max)
void Print (const char *filename, int options=0, double ovalue=0., int ww=0, int wh=0)
Member functions (get)



TH1D * GetHistogram ()
double GetIntegral (double valuemin, double valuemax)
double GetLimit (double probability)
double GetMean ()
double GetMedian ()
double GetMode ()
double GetPValue (double probability)
double GetQuantile (double probablity)
double GetRMS ()
Member functions (set)



void SetDefaultCLLimit (double limit)
void SetGlobalMode (double mode)
void SetHistogram (TH1D *hist)

Private Attributes

double fDefaultCLLimit
TH1D * fHistogram
double fMode
int fModeFlag

Detailed Description

A class for handling 1D distributions.

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

Definition at line 30 of file BCH1D.h.


Constructor & Destructor Documentation

BCH1D::BCH1D (  ) 

The default constructor.

Definition at line 31 of file BCH1D.cxx.

00032 {
00033    fHistogram = 0;
00034    fDefaultCLLimit = 95.; // in percent
00035 
00036    fModeFlag = 0;
00037 }

BCH1D::BCH1D ( TH1D *  hist  )  [inline]

The default constructor.

Definition at line 43 of file BCH1D.h.

00044          { fHistogram = hist; };

BCH1D::~BCH1D (  ) 

The default destructor.

Definition at line 41 of file BCH1D.cxx.

00042 {
00043    if (fHistogram) delete fHistogram;
00044 }


Member Function Documentation

void BCH1D::Draw ( int  options = 0,
double  ovalue = 0. 
)

Draw distribution into the active canvas.

Parameters:
options Drawing options:
0 = band mode [default],
1 = draw vertical line,
2 = band mode with minimal interval
ovalue Option specific value. For option 0, if ovalue is nonzero a limit is to be drawn rather than central band with ovalue being the per cent value of the limit. If negative, limit is drawn from minimum, if positive limit is drawn from maximum. Allowed values are 68 < |limit| < 100. If mode is outside the band, the limit is drawn automatically. The default limit can be changed by BCH1D::SetDefaultCLLimit(int limit).
For option 1 the ovalue defines where the line is drawn.
For option 2 the ovalue sets the content of the minimal interval in per cent. If omitted a 68% minimal interval will be drawn.

Definition at line 137 of file BCH1D.cxx.

00138 {
00139    double min, max;
00140    double mode;
00141    double thismode = this->GetMode();
00142 
00143    int nbins = fHistogram->GetNbinsX();
00144 
00145    fHistogram->Scale(1./fHistogram->Integral("width"));
00146 
00147    if(fModeFlag)
00148       mode=fMode;
00149    else
00150       mode=thismode;
00151 
00152    // define temporary line.
00153    TLine * line;
00154 
00155    // check drawing option.
00156    switch(options)
00157    {
00158       // Draw a band between 16% and 84% probability.
00159       // If the mode is outside the band only draw a limit.
00160       case 0:
00161 
00162          if (fabs(ovalue) >= 100 || fabs(ovalue) < 68)
00163          {
00164             min = this -> GetQuantile(.16);
00165             max = this -> GetQuantile(.84);
00166 
00167 
00168 //          if ( thismode > max || fHistogram -> FindBin(thismode) == fHistogram -> GetNbinsX() )
00169             if ( fHistogram -> FindBin(thismode) == fHistogram -> GetNbinsX() )
00170             {
00171                min = this -> GetQuantile(1.-(double)fDefaultCLLimit/100.);
00172                max = fHistogram->GetXaxis()->GetXmax();
00173                ovalue = fDefaultCLLimit;
00174             }
00175 //          else if (thismode < min || fHistogram->FindBin(thismode)==1)
00176             else if ( fHistogram->FindBin(thismode)==1)
00177             {
00178                min = fHistogram->GetXaxis()->GetXmin();
00179                max = this -> GetQuantile((double)fDefaultCLLimit/100.);
00180                ovalue = -fDefaultCLLimit;
00181             }
00182          }
00183 
00184          else if(ovalue < 0)
00185          {
00186             min = fHistogram->GetXaxis()->GetXmin();
00187             max = this -> GetQuantile(-ovalue/100.);
00188          }
00189          else
00190          {
00191             min = this -> GetQuantile((1-ovalue)/100.);
00192             max = fHistogram->GetXaxis()->GetXmax();
00193          }
00194 
00195          // do the drawing
00196          this->DrawShadedLimits(mode, min, max, ovalue);
00197 
00198          break;
00199 
00200       // Draw a line at "ovalue".
00201       case 1:
00202 
00203          fHistogram -> Draw();
00204          min = fHistogram->GetBinLowEdge(1);
00205          max = fHistogram->GetBinLowEdge(nbins+1);
00206          if(min<=ovalue && ovalue<=max)
00207          {
00208             line = new TLine();
00209             line -> SetLineColor(kRed);
00210             line -> DrawLine(ovalue, 0., ovalue, 1.05 * fHistogram -> GetMaximum());
00211          }
00212 
00213          break;
00214 
00215       // Draw a shaded band at the smallest interval.
00216       case 2:
00217 
00218          if(ovalue<50) // just to ensure there's some sense in the number
00219             ovalue = 68.; // default is 68%
00220 
00221          this->GetSmallestInterval(min, max, ovalue/100.);
00222          this->DrawShadedLimits(mode, min, max, 0.);
00223 
00224          break;
00225 
00226       // Draw a shaded band at the smallest intervals
00227       case 3:
00228 
00229          if(ovalue<50) // just to ensure there's some sense in the number
00230             ovalue = 68.; // default is 68%
00231 
00232          this -> DrawSmallest(mode,ovalue);
00233 
00234          break;
00235 
00236       // Sort out bad options and warn.
00237       default:
00238 
00239          BCLog::Out(BCLog::warning, BCLog::warning, Form("BCH1D::Draw : Invalid option %d",options));
00240          break;
00241    }
00242 }

void BCH1D::DrawShadedLimits ( double  mode,
double  min,
double  max,
double  limit = 0 
)

Draw distribution with band between min and max and with marker at the mode. Write the location of the mode with uncertainties. If limit is specified, draw CL limit. Allowed values are 68 < |limit| < 100.

Definition at line 246 of file BCH1D.cxx.

00247 {
00248    double maximum = fHistogram -> GetMaximum();
00249 
00250    double x0 = mode;
00251    double y0 = fHistogram->GetBinContent( fHistogram->FindBin(mode) );
00252 
00253    double x1 = fHistogram->GetMean();
00254    double y1 = fHistogram->GetBinContent( fHistogram->FindBin(x1) );
00255 
00256    double x2 = this -> GetQuantile(.5); // median
00257    double y2 = fHistogram->GetBinContent( fHistogram->FindBin(x2) );
00258 
00259    double ysize = maximum*1.3;
00260 
00261    double xmin = fHistogram->GetXaxis()->GetXmin();
00262    double xmax = fHistogram->GetXaxis()->GetXmax();
00263 
00264    // draw histogram with axes first
00265    TH2D * hsc = new TH2D(
00266          TString::Format("h2scale_%s_%d",fHistogram->GetName(),BCLog::GetHIndex()),"",
00267          10, xmin, xmax, 10, 0., ysize);
00268    hsc -> SetStats(0);
00269    hsc -> SetXTitle(fHistogram->GetXaxis()->GetTitle());
00270    hsc -> SetYTitle(fHistogram->GetYaxis()->GetTitle());
00271    hsc -> Draw();
00272 
00273    // draw histogram
00274    fHistogram -> SetLineWidth(1);
00275    fHistogram -> Draw("same");
00276 
00277    // draw yellow shaded region between min and max
00278    TH1D * hist_shaded = this->GetSubHisto(min,max,TString::Format("%s_sub_%d",fHistogram->GetName(),BCLog::GetHIndex()));
00279    hist_shaded -> SetFillStyle(1001);
00280    hist_shaded -> SetFillColor(kYellow);
00281 
00282    // draw shaded histogram
00283    hist_shaded -> Draw("same");
00284 
00285    gPad->RedrawAxis();
00286 
00287    // draw triangle for mode
00288    TPolyLine * tmax;
00289 
00290    double dx = 0.01*(xmax-xmin);
00291    double dy = 0.04*(ysize);
00292    y0+=dy/5.;
00293    double tmax_x[] = {x0, x0 + dx, x0 - dx, x0};
00294    double tmax_y[] = {y0, y0 + dy, y0 + dy, y0};
00295    tmax = new TPolyLine(4,tmax_x,tmax_y);
00296    tmax->SetLineColor(kRed);
00297    tmax->SetLineWidth(1);
00298    tmax->SetFillColor(kRed);
00299    tmax->Draw();
00300    tmax->Draw("f");
00301 
00302    // draw triangle for mean
00303    // draw diamond for mean
00304    TPolyLine * tmean;
00305 
00306    y1+=dy/5.;
00307 // double tmean_x[] = {x1, x1 + dx, x1 - dx, x1};
00308 // double tmean_y[] = {y1, y1 + dy, y1 + dy, y1};
00309    double tmean_x[] = {x1, x1 + dx, x1 , x1 - dx, x1};
00310    double tmean_y[] = {y1, y1 + dy/2., y1 + dy, y1 + dy/2., y1};
00311    tmean = new TPolyLine(5,tmean_x,tmean_y);
00312    tmean->SetLineColor(kBlue);
00313 // tmean->SetFillColor(kWhite);
00314    tmean->SetLineWidth(1);
00315 // tmean->SetLineStyle(1);
00316    tmean->Draw();
00317 
00318 /*
00319    // draw arrow for median
00320    TPolyLine * tmed;
00321    TPolyLine * tmed2;
00322 
00323 // y2+=dy/5.;
00324    y2=0.+dy/5.;
00325    double tmed_x[] = {x2 + dx, x2, x2 - dx};
00326    double tmed_y[] = {y2 + dy, y2, y2 + dy};
00327    double tmed2_x[] = {x2, x2};
00328    double tmed2_y[] = {y2, y2 + dy*2.};
00329    tmed = new TPolyLine(3,tmed_x,tmed_y);
00330    tmed2 = new TPolyLine(2,tmed2_x,tmed2_y);
00331    tmed->SetLineColor(kGreen+2);
00332    tmed->SetLineWidth(1);
00333    tmed2->SetLineColor(kGreen+2);
00334    tmed2->SetLineWidth(1);
00335    tmed->Draw();
00336    tmed2->Draw();
00337 */
00338    // draw dashed line for median
00339    TLine * line;
00340    line = new TLine();
00341    line -> SetLineStyle(2);
00342    line -> SetLineColor(kGreen+2);
00343    line -> DrawLine(x2, 0., x2, y2);
00344 
00345    // write mode location and shaded band
00346 
00347    // format of the number
00348    double delta_max = fmax(fabs(max-x1),fabs(x1-min));
00349 
00350    int sd = 2 + (int)log10(fabs(x1/delta_max));
00351 
00352    if( (int)log10(x1) > (int)log10(delta_max) )
00353       sd++;
00354 
00355    TLatex * tmax_text = new TLatex();
00356    tmax_text->SetTextSize(0.045);
00357    tmax_text->SetTextFont(62);
00358    tmax_text->SetTextAlign(22); // center
00359 // tmax_text->SetTextAlign(13); // top-left
00360 
00361    double xprint=(xmax+xmin)/2.;
00362 // double xprint=xmin+(xmax-xmin)/20.;
00363    double yprint=ysize*(1-1.4*tmax_text->GetTextSize());
00364 
00365    if(fabs(limit)<50) // just to ensure there's some sense in the number
00366       tmax_text->DrawLatex(xprint,yprint,
00367          TString::Format( TString::Format("%%s^{med} = %%.%dg ^{+%%.2g}_{ -%%.2g}",sd),
00368             fHistogram->GetXaxis()->GetTitle(), x2, max-x2, x2-min));
00369 
00370    else if (limit<0)
00371       tmax_text->DrawLatex(xprint,yprint,
00372          TString::Format( TString::Format("%%s (%d%%%% prob.) < %%.4g",-(int)limit),
00373             fHistogram->GetXaxis()->GetTitle(), max));
00374 
00375    else if (limit>0)
00376       tmax_text->DrawLatex(xprint,yprint,
00377          TString::Format( TString::Format("%%s (%d%%%% prob.) > %%.4g",(int)limit),
00378             fHistogram->GetXaxis()->GetTitle(), min));
00379 
00380 /*
00381    TLegend * leg = new TLegend(.61,.7,.9,.88);
00382    leg -> SetFillColor(kWhite);
00383    leg -> SetFillStyle(0);
00384    leg -> SetBorderSize(0);
00385    leg -> AddEntry(line,"Median", "l");
00386    TH1D * hh0 = new TH1D("hh0","",10,0,10);
00387    hh0->SetMarkerStyle(23);
00388    hh0->SetMarkerColor(kBlue);
00389    hh0->SetMarkerSize(1);
00390    leg -> AddEntry(hh0,"Mean","p");
00391    TH1D * hh1 = new TH1D("hh1","",10,0,10);
00392    hh1->SetMarkerStyle(23);
00393    hh1->SetMarkerColor(kRed);
00394    hh1->SetMarkerSize(1);
00395    leg -> AddEntry(hh1,"Mode","p");
00396    leg -> AddEntry(hist_shaded, "Central 68%", "f");
00397    leg ->Draw();
00398 */
00399 
00400 }

void BCH1D::DrawSmallest ( double  mode,
double  prob,
bool  drawmean = true 
)

Draw distribution with bands so that the total shaded area is the smallest possible containing and integral of "prob". Draw the location of the mean and median if requested (default).

Definition at line 404 of file BCH1D.cxx.

00405 {
00406    // prepare triangle for mode
00407    TPolyLine * tmax;
00408 
00409    double x0 = mode;
00410    double y0 = fHistogram->GetBinContent( fHistogram->FindBin(mode) );
00411    double xmin = fHistogram->GetXaxis()->GetXmin();
00412    double xmax = fHistogram->GetXaxis()->GetXmax();
00413    double ysize = 1.3 * fHistogram -> GetMaximum();
00414 
00415    double x1 = fHistogram->GetMean();
00416    double y1 = fHistogram->GetBinContent( fHistogram->FindBin(x1) );
00417 
00418    double x2 = this -> GetQuantile(.5); // median
00419    double y2 = fHistogram->GetBinContent( fHistogram->FindBin(x2) );
00420 
00421    double dx = 0.01*(xmax-xmin);
00422    double dy = 0.04*(ysize);
00423    double tmax_x[] = {x0, x0 + dx, x0 - dx, x0};
00424    double tmax_y[] = {y0, y0 + dy, y0 + dy, y0};
00425    tmax = new TPolyLine(4,tmax_x,tmax_y);
00426    tmax->SetLineColor(kRed);
00427    tmax->SetFillColor(kRed);
00428 
00429    // draw histogram with axes first
00430    TH2D * hsc = new TH2D(
00431          TString::Format("h2scale_%s_%d",fHistogram->GetName(),BCLog::GetHIndex()),"",
00432          10, xmin, xmax, 10, 0., ysize);
00433    hsc -> SetStats(0);
00434    hsc -> SetXTitle(fHistogram->GetXaxis()->GetTitle());
00435    hsc -> SetYTitle(fHistogram->GetYaxis()->GetTitle());
00436    hsc -> Draw();
00437 
00438    // histogram to be filled with band
00439    TH1D * hist_temp1 = (TH1D*) fHistogram -> Clone();
00440    hist_temp1 -> Scale(1.0/fHistogram -> Integral("width"));
00441    hist_temp1 -> SetFillColor(kYellow);
00442    hist_temp1 -> SetFillStyle(1001);
00443 
00444    // temporary histogram
00445    TH1D * hist_temp2 = (TH1D*) fHistogram -> Clone();
00446    hist_temp2 -> Scale(1.0/fHistogram -> Integral("width"));
00447 
00448    // clear content
00449    hist_temp1 -> Reset();
00450 
00451    // loop over original histogram and copy bin untils integral equals
00452    // "prob"
00453    prob /= 100.;
00454 
00455    // integral
00456    double sum = 0.0;
00457 
00458    int n=0;
00459    int nbins=fHistogram->GetNbinsX();
00460 
00461    // loop
00462    while (sum < prob && n < nbins)
00463    {
00464       n++;
00465 
00466       // find bin with maximum
00467       int bin = hist_temp2 -> GetMaximumBin();
00468 
00469       // copy bin to new histogram
00470       double val = hist_temp2 -> GetBinContent(bin);
00471       hist_temp1 -> SetBinContent(bin, val);
00472 
00473       // remove maximum from temporary histogram
00474       hist_temp2 -> SetBinContent(bin, 0.0);
00475 
00476       // integrate by summing
00477       sum += val * hist_temp2->GetBinWidth(bin);
00478    }
00479 
00480    // scale histogram
00481    hist_temp1 -> Scale(fHistogram -> Integral("width"));
00482 
00483    // draw histograms
00484    fHistogram -> Draw("same");
00485    hist_temp1 -> Draw("same");
00486 
00487    // draw triangle for mode
00488    tmax->Draw("f");
00489 
00490    if(drawmean)
00491    {
00492       // draw triangle for mean
00493       // draw diamond for mean
00494       TPolyLine * tmean;
00495 
00496       y1+=dy/5.;
00497 //    double tmean_x[] = {x1, x1 + dx, x1 - dx, x1};
00498 //    double tmean_y[] = {y1, y1 + dy, y1 + dy, y1};
00499       double tmean_x[] = {x1, x1 + dx, x1 , x1 - dx, x1};
00500       double tmean_y[] = {y1, y1 + dy/2., y1 + dy, y1 + dy/2., y1};
00501       tmean = new TPolyLine(5,tmean_x,tmean_y);
00502       tmean->SetLineColor(kBlue);
00503 //    tmean->SetFillColor(kWhite);
00504       tmean->SetLineWidth(1);
00505 //    tmean->SetLineStyle(1);
00506       tmean->Draw();
00507 
00508       // draw dashed line for median
00509       TLine * line;
00510       line = new TLine();
00511       line -> SetLineStyle(2);
00512       line -> SetLineColor(kGreen+2);
00513       line -> DrawLine(x2, 0., x2, y2);
00514    }
00515 
00516    // free memory
00517    delete hist_temp2;
00518 }

TH1D* BCH1D::GetHistogram (  )  [inline]
Returns:
The one-dimensional histogram.

Definition at line 57 of file BCH1D.h.

00058          { return fHistogram; };

double BCH1D::GetIntegral ( double  valuemin,
double  valuemax 
)

Returns the integral of distribution the between two values.

Parameters:
valuemin The value from which the intergration is done.
valuemax The value up to which the intergration is done.
Returns:
The integral.

Definition at line 71 of file BCH1D.cxx.

00072 {
00073    double integral = 0;
00074 
00075    int binmin = fHistogram -> FindBin(valuemin);
00076    int binmax = fHistogram -> FindBin(valuemax);
00077 
00078    // use ROOT function to calculate integral.
00079    integral = fHistogram -> Integral(binmin, binmax);
00080 
00081    return integral;
00082 }

double BCH1D::GetLimit ( double  probability  )  [inline]

Return the quantily of the distribution

Parameters:
probability The probability.
Returns:
The quantile of the distribution for the probability.
See also:
GetQuantile(double probablity)

Definition at line 86 of file BCH1D.h.

00087          { return this -> GetQuantile(probability); };

double BCH1D::GetMean (  )  [inline]
Returns:
The mean of the distribution.

Definition at line 62 of file BCH1D.h.

00063          { return fHistogram -> GetMean(); };

double BCH1D::GetMedian (  )  [inline]
Returns:
The median of the distribution.

Definition at line 71 of file BCH1D.h.

00072          { return this -> GetQuantile(0.5); };

double BCH1D::GetMode (  ) 
Returns:
The mode of the distribution.

Definition at line 48 of file BCH1D.cxx.

00049 {
00050    return fHistogram -> GetBinCenter(fHistogram -> GetMaximumBin());
00051 }

double BCH1D::GetPValue ( double  probability  ) 

Returns the p-value. Returns the integral from 0 to the probability.

Parameters:
probability Upper limit of integration.
Returns:
The p-value.

Definition at line 86 of file BCH1D.cxx.

00087 {
00088    // use ROOT function to calculate the integral from 0 to "probability".
00089    return fHistogram -> Integral(1, fHistogram -> FindBin(probability));
00090 }

double BCH1D::GetQuantile ( double  probablity  ) 

Returns the quantile of the distribution.

Parameters:
probability The probability.
Returns:
The quantile of the distribution for the probability.
See also:
GetLimit(double probability)

Definition at line 55 of file BCH1D.cxx.

00056 {
00057    int nquantiles = 1;
00058    double quantiles[1];
00059    double probsum[1];
00060 
00061    probsum[0] = probability;
00062 
00063    // use ROOT function to calculat quantile.
00064    fHistogram -> GetQuantiles(nquantiles, quantiles, probsum);
00065 
00066    return quantiles[0];
00067 }

double BCH1D::GetRMS (  )  [inline]
Returns:
The RMS of the distribution.

Definition at line 91 of file BCH1D.h.

00092          { return fHistogram -> GetRMS(); };

double BCH1D::GetSmallestInterval ( double &  min,
double &  max,
double  content = .68 
)

Calculate the minimal interval of the distribution containing a given content.

Parameters:
min calculated minimum of the interval
max calculated maximum of the interval
content content of the interval [default is .68]
Returns:
the content of the histogram between min and max

Definition at line 522 of file BCH1D.cxx.

00523 {
00524    if(content<=0. || content>= 1.)
00525       return 0.;
00526 
00527    int nbins = fHistogram -> GetNbinsX();
00528 
00529    double factor = fHistogram->Integral("width");
00530    if(factor==0)
00531       return 0.;
00532 
00533    // normalize if not yet done
00534    fHistogram->Scale(1./factor);
00535 
00536    double xmin = fHistogram->GetXaxis()->GetXmin();
00537    double xmax = fHistogram->GetXaxis()->GetXmax();
00538    double width = xmax - xmin;
00539 
00540    double xdown=xmin;
00541    double xup=xmax;
00542 
00543    int ndiv = 100;
00544    if(nbins<100)
00545       ndiv = 1000;
00546    if(nbins>1000)
00547       ndiv = 10;
00548 
00549    int warn=0;
00550 
00551    double integral_out=0.;
00552 
00553    // loop through the bins
00554    for(int i=1;i<nbins+1;i++)
00555    {
00556       if(fHistogram->Integral(i,nbins,"width") < content)
00557          break;
00558 
00559       // get width of the starting bin
00560       double firstbinwidth = fHistogram->GetBinWidth(i);
00561 
00562       // divide i-th bin in ndiv divisions
00563       // integrate starting at each of these divisions
00564       for(int j=0;j<ndiv;j++)
00565       {
00566          double dxdown = (double)(ndiv-j)/(double)ndiv * firstbinwidth;
00567          xdown = fHistogram->GetBinLowEdge(i) + firstbinwidth - dxdown;
00568          double integral = dxdown * fHistogram->GetBinContent(i);
00569 
00570          if(integral>content)
00571          {
00572             // content fits within one bin
00573             xup = xdown + content / fHistogram->GetBinContent(i);
00574             warn = 1;
00575          }
00576          else
00577             for(int k=i+1;k<nbins+1;k++)
00578             {
00579 
00580                double thisbin = fHistogram->GetBinContent(k) * fHistogram->GetBinWidth(k);
00581 
00582                if(integral+thisbin==content)
00583                {
00584                   xup = fHistogram->GetBinLowEdge(k+1);
00585                   break;
00586                }
00587 
00588                if(integral+thisbin>content)
00589                {
00590                   xup = fHistogram->GetBinLowEdge(k) + (content-integral)/thisbin * fHistogram->GetBinWidth(k);
00591                   integral += thisbin * (xup-fHistogram->GetBinLowEdge(k))/fHistogram->GetBinWidth(k);
00592                   break;
00593                }
00594 
00595                integral += thisbin;
00596             }
00597 
00598          if(integral < content)
00599             continue;
00600 
00601          if(xup - xdown < width)
00602          {
00603             // store necessary information
00604             width = xup - xdown;
00605             xmin  = xdown;
00606             xmax  = xup;
00607             integral_out = integral;
00608          }
00609       }
00610    }
00611 
00612    if(warn)
00613    {
00614       BCLog::Out(BCLog::warning,BCLog::warning,
00615             Form("BCH1D::GetSmallestInterval : The requested content of %d%% fits within one bin.",(int)(content*100)));
00616       BCLog::Out(BCLog::warning,BCLog::warning,
00617             "BCH1D::GetSmallestInterval : MAKE FINER BINNING (OR CHANGE RANGE) !!!");
00618    }
00619 
00620    // restore normalization to state before calling this method
00621    fHistogram->Scale(factor);
00622 
00623    min=xmin;
00624    max=xmax;
00625 
00626    return integral_out;
00627 }

TH1D * BCH1D::GetSmallestIntervalHistogram ( double  level  ) 

Definition at line 749 of file BCH1D.cxx.

00750 {
00751    // create a new histogram which will be returned and all yellow
00752    TH1D * hist_yellow = (TH1D*) fHistogram -> Clone();
00753    hist_yellow -> Reset();
00754    hist_yellow -> SetFillStyle(1001);
00755    hist_yellow -> SetFillColor(kYellow);
00756 
00757    // copy a temporary histogram
00758    TH1D * hist_temp = (TH1D*) fHistogram -> Clone(TString::Format("%s_%i",fHistogram->GetName(),BCLog::GetHIndex()));
00759    double factor = hist_temp -> Integral("");
00760 
00761    if(factor == 0)
00762       return 0;
00763 
00764    hist_temp -> Scale(1. / factor);
00765 
00766    // here's the algorithm:
00767    // 1. find the maximum bin in the temporary histogram and copy it to
00768    // the yellow histogram.
00769    // 2. remove this bin from the temporary histogram.
00770    // 3. repeat this until a total of "level" probability is copied to
00771    // the yellow histogram.
00772 
00773    double sumprob = 0.0;
00774 
00775    while (sumprob < level)
00776    {
00777       // find maximum bin and value
00778       int bin = hist_temp -> GetMaximumBin();
00779       double value = hist_temp -> GetMaximum();
00780 
00781       // copy "1" into the corresponding bin in the yellow histogram
00782       hist_yellow -> SetBinContent(bin, 1.0);
00783 
00784       // set the bin value in the temporary histogram to zero
00785       hist_temp -> SetBinContent(bin, 0.0);
00786 
00787       // increase probability sum
00788       sumprob += value;
00789    }
00790 
00791    delete hist_temp;
00792 
00793    return hist_yellow;
00794 }

std::vector< double > BCH1D::GetSmallestIntervals ( double  content = 0.68  ) 

Definition at line 798 of file BCH1D.cxx.

00799 {
00800    std::vector <double> v;
00801 
00802    TH1D * hist = this -> GetSmallestIntervalHistogram(content);
00803 
00804    int nbins = hist -> GetNbinsX();
00805    int ninter = 0;
00806    int lastbin = -1;
00807 
00808    double max = -1;
00809    double localmax = -1;
00810    double localmaxpos = -1;
00811    double localint = 0;
00812    bool flag = false;
00813 
00814    for (int i = 1; i <= nbins; ++i)
00815    {
00816       // interval starts here
00817       if (!flag && hist -> GetBinContent(i) > 0.)
00818       {
00819          flag = true;
00820          v.push_back(hist -> GetBinLowEdge(i));
00821 
00822          // remember start position of the interval
00823          lastbin = i;
00824 
00825          // increase number of intervals
00826          ninter++;
00827 
00828          // reset local maximum
00829          localmax = fHistogram -> GetBinContent(i);
00830          localmaxpos = hist -> GetBinLowEdge(i);
00831 
00832          // reset local integral
00833          localint = 0;
00834       }
00835 
00836       // interval stops here
00837       if ((flag && !(hist -> GetBinContent(i) > 0.)) || (flag && i == nbins))
00838       {
00839          flag = false;
00840          v.push_back(hist -> GetBinLowEdge(i) + hist -> GetBinWidth(i));
00841 
00842          // set right bin to maximum if on edge
00843          if (i == nbins && localmax < fHistogram -> GetBinContent(i))
00844             localmaxpos = hist -> GetBinCenter(i) + 0.5 * hist -> GetBinWidth(i);
00845 
00846          // find the absolute maximum
00847          if (localmax > max)
00848             max = localmax;
00849 
00850          // save local maximum
00851          v.push_back(localmax);
00852          v.push_back(localmaxpos);
00853 
00854          // save local integral
00855          v.push_back(localint);
00856       }
00857 
00858       // find local maximum
00859       if (i < nbins && localmax < fHistogram -> GetBinContent(i))
00860       {
00861          localmax = fHistogram -> GetBinContent(i);
00862          localmaxpos = hist -> GetBinCenter(i);
00863       }
00864 
00865       // increase area
00866       localint += fHistogram -> GetBinContent(i) / fHistogram -> Integral();
00867    }
00868 
00869    // rescale absolute heights to relative heights
00870    for (int i = 0; i < ninter; ++i)
00871       v[i*5+2] = v.at(i*5+2) / max;
00872 
00873    return v;
00874 }

TH1D * BCH1D::GetSubHisto ( double  min,
double  max,
const char *  name 
)

Get histogram with bins outside min, max band being zero. The new histogram can have 2 more bins than the original one as the bins where min and max fall into will be split in two (except for the case when min and/or max are equal to some of the original bin boundaries.

Parameters:
min lower boundary of the non-zero interval
max upper boundary of the non-zero interval
Returns:
new histogram which is nonzero only between min and max

Definition at line 670 of file BCH1D.cxx.

00671 {
00672    if(min>max)
00673    {
00674       double t=min;
00675       min=max;
00676       max=t;
00677    }
00678 
00679    int imin = fHistogram->FindBin(min);
00680    int imax = fHistogram->FindBin(max);
00681 
00682    int nbins = fHistogram->GetNbinsX();
00683    double xmin = fHistogram->GetXaxis()->GetXmin();
00684    double xmax = fHistogram->GetXaxis()->GetXmax();
00685 
00686    if( min==max || (min<=xmin && max>=xmax) )
00687    {
00688       TH1D * h0 = (TH1D*) fHistogram->Clone(name);
00689       return h0;
00690    }
00691 
00692    if (imin<1)
00693    {
00694       imin=1;
00695       min=xmin;
00696    }
00697    if (imax>nbins)
00698    {
00699       imax=nbins;
00700       max=xmax;
00701    }
00702 
00703    double * xb = new double[nbins+3]; // nbins+1 original bin edges + 2 new bins
00704    int n=0; // counter
00705 
00706    int domin=1;
00707    int domax=1;
00708 
00709    for (int i=1;i<nbins+2;i++)
00710    {
00711       double x0 = fHistogram->GetBinLowEdge(i);
00712 
00713       if (min<x0 && domin)
00714       {
00715          xb[n++]=min;
00716          domin=0;
00717       }
00718       else if (min==x0)
00719          domin=0;
00720 
00721       if (max<x0 && domax)
00722       {
00723          xb[n++]=max;
00724          domax=0;
00725       }
00726       else if (max==x0)
00727          domax=0;
00728 
00729       xb[n++]=x0;
00730    }
00731 
00732    // now define the new histogram
00733    TH1D * h0 = new TH1D(name,"",n-1,xb);
00734    for(int i=1;i<n;i++)
00735    {
00736       double x0 = h0->GetBinCenter(i);
00737       if(x0<min || x0>max)
00738          continue;
00739 
00740       int bin=fHistogram->FindBin(x0);
00741       h0->SetBinContent(i, fHistogram->GetBinContent(bin));
00742    }
00743 
00744    return h0;
00745 }

double BCH1D::IntegralWidth ( double  min,
double  max 
)

Calculate integral of the distribution between min and max.

Parameters:
min lower boundary of the integrated interval
max upper boundary of the integrated interval
Returns:
integral calculated as sum of BinContent*BinWidth

Definition at line 631 of file BCH1D.cxx.

00632 {
00633    int imin = fHistogram->FindBin(min);
00634    int imax = fHistogram->FindBin(max);
00635 
00636    int nbins = fHistogram->GetNbinsX();
00637 
00638    // if outside of histogram range, return -1.
00639    if ( imin<1 || imin>nbins || imax<1 || imax>nbins )
00640       return -1.;
00641 
00642    if ( imin==imax )
00643       return -1.;
00644 
00645    // swap values if necessary
00646    if (imin>imax)
00647    {
00648       int i=imin;
00649       double x=min;
00650       imin=imax, min=max;
00651       imax=i, max=x;
00652    }
00653 
00654    // calculate first bin
00655    double first = ( fHistogram->GetBinLowEdge(imin+1) - min ) * fHistogram->GetBinContent(imin);
00656 
00657    // calculate last bin
00658    double last = ( max - fHistogram->GetBinLowEdge(imax) ) * fHistogram->GetBinContent(imax);
00659 
00660    // calculate inbetween
00661    double inbetween=0.;
00662    if(imax-imin>1)
00663       inbetween = fHistogram->Integral(imin+1, imax-1, "width");
00664 
00665    return first + inbetween + last;
00666 }

void BCH1D::Print ( const char *  filename,
int  options = 0,
double  ovalue = 0.,
int  ww = 0,
int  wh = 0 
)

Print distribution into a PostScript file.

Parameters:
filename Output 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 parameters options and ovalue look at BCH1D::Draw() method.

Definition at line 108 of file BCH1D.cxx.

00109 {
00110    char file[256];
00111    int i=0;
00112    while(i<255 && *filename!='\0')
00113       file[i++]=*filename++;
00114    file[i]='\0';
00115 
00116 // temporary
00117    fHistogram -> SetLineWidth(1);
00118 
00119    // create temporary canvas
00120    TCanvas * c;
00121    if(ww > 0 && wh > 0)
00122       c = new TCanvas("c","c",ww,wh);
00123    else
00124       c = new TCanvas("c","c");
00125 
00126    c -> cd();
00127    this->Draw(options, ovalue);
00128 
00129    gPad->RedrawAxis();
00130 
00131    // print to file.
00132    c -> Print(file);
00133 }

void BCH1D::SetDefaultCLLimit ( double  limit  ) 

Set default probability limits. Allowed values are between 68% and 100%. The default value is 95%.

Definition at line 94 of file BCH1D.cxx.

00095 {
00096    // check if limit is between 68% and 100%. Give a warning if not ...
00097    if(limit>=100. || limit<68.)
00098       BCLog::Out(BCLog::warning,BCLog::warning,
00099          Form("BCH1D::SetDefaultCLLimit : Value %.1f out of range. Keeping default %.1f%% CL limit.",limit,fDefaultCLLimit));
00100 
00101    // ... or set value if true.
00102    else
00103       fDefaultCLLimit = limit;
00104 }

void BCH1D::SetGlobalMode ( double  mode  )  [inline]

Set global mode

Definition at line 125 of file BCH1D.h.

00126          { fMode=mode; fModeFlag=1; };

void BCH1D::SetHistogram ( TH1D *  hist  )  [inline]

Sets the histogram.

Definition at line 115 of file BCH1D.h.

00116          { fHistogram = hist; };


Member Data Documentation

double BCH1D::fDefaultCLLimit [private]

Default confidence level limit

Definition at line 214 of file BCH1D.h.

TH1D* BCH1D::fHistogram [private]

The 1D histogram

Definition at line 210 of file BCH1D.h.

double BCH1D::fMode [private]

Global mode

Definition at line 218 of file BCH1D.h.

int BCH1D::fModeFlag [private]

"Is there a global mode?" flag

Definition at line 222 of file BCH1D.h.


The documentation for this class was generated from the following files:

Generated on Tue Oct 6 09:48:21 2009 for Bayesian Analysis Toolkit by  doxygen 1.6.1