Private Attributes

BCH1D Class Reference

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

#include <BCH1D.h>

List of all members.

Public Member Functions

Constructors and destructors

 BCH1D ()
 BCH1D (TH1D *hist)
 ~BCH1D ()
Member functions (get)

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

void SetHistogram (TH1D *hist)
void SetDefaultCLLimit (double limit)
void SetGlobalMode (double mode)
Member functions (miscellaneous methods)

void Print (const char *filename, int options=0, double ovalue=0., int ww=0, int wh=0)
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)
void DrawLegend (const char *text)
double GetSmallestInterval (double &min, double &max, double content=.68)
TH1D * GetSmallestIntervalHistogram (double level)
std::vector< double > GetSmallestIntervals (double content=0.68)
double IntegralWidth (double min, double max)
TH1D * GetSubHisto (double min, double max, const char *name)

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.

{
   fHistogram = 0;
   fDefaultCLLimit = 95.; // in percent

   fModeFlag = 0;
}

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

The default constructor.

Definition at line 43 of file BCH1D.h.

         { fHistogram = hist; };

BCH1D::~BCH1D (  ) 

The default destructor.

Definition at line 41 of file BCH1D.cxx.

{
   if (fHistogram) delete fHistogram;
}


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.

{
   double min, max;
   double mode;
   double thismode = this->GetMode();

   int nbins = fHistogram->GetNbinsX();

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

   if(fModeFlag)
      mode=fMode;
   else
      mode=thismode;

   // define temporary line.
   TLine * line;

   //control if legend is required
   bool flagLegend=false;
   char confidenceLabel[25];

   // check drawing option.
   switch(options)
   {
      // Draw a band between 16% and 84% probability.
      // If the mode is outside the band only draw a limit.
      case 0:
         if (fabs(ovalue) >= 100 || ovalue==0.)
         {//default case if no args to Draw() supplied

            min = this -> GetQuantile(.16);
            max = this -> GetQuantile(.84);

            //draw a legend later
            flagLegend = true;
            sprintf(confidenceLabel, "Central 68%%");

//          if ( thismode > max || fHistogram -> FindBin(thismode) == fHistogram -> GetNbinsX() )
            if ( fHistogram -> FindBin(thismode) == fHistogram -> GetNbinsX() )
            {
               min = this -> GetQuantile(1.-(double)fDefaultCLLimit/100.);
               max = fHistogram->GetXaxis()->GetXmax();
               ovalue = fDefaultCLLimit;
               sprintf(confidenceLabel, "%g%% region", fDefaultCLLimit);
            }
//          else if (thismode < min || fHistogram->FindBin(thismode)==1)
            else if ( fHistogram->FindBin(thismode)==1)
            {
               min = fHistogram->GetXaxis()->GetXmin();
               max = this -> GetQuantile((double)fDefaultCLLimit/100.);
               ovalue = -fDefaultCLLimit;
               sprintf(confidenceLabel, "%g%% region", fDefaultCLLimit);
            }
         }

         else if(ovalue < 0)
         {
            min = fHistogram->GetXaxis()->GetXmin();
            max = this -> GetQuantile(-ovalue/100.);
         }
         else
         {
            min = this -> GetQuantile(1. - ovalue/100.);
            max = fHistogram->GetXaxis()->GetXmax();
         }

         // do the drawing
         this->DrawShadedLimits(mode, min, max, ovalue);

         // add legend for the symbols mean, mode, median, confidence band
         if(flagLegend)
            this ->DrawLegend(confidenceLabel);


         break;

      // Draw a line at "ovalue".
      case 1:

         fHistogram -> Draw();
         min = fHistogram->GetBinLowEdge(1);
         max = fHistogram->GetBinLowEdge(nbins+1);
         if(min<=ovalue && ovalue<=max)
         {
            line = new TLine();
            line -> SetLineColor(kRed);
            line -> DrawLine(ovalue, 0., ovalue, 1.05 * fHistogram -> GetMaximum());
         }

         break;

      // Draw a shaded band at the smallest interval.
      case 2:

         if(ovalue<50) // just to ensure there's some sense in the number
            ovalue = 68.; // default is 68%

         this->GetSmallestInterval(min, max, ovalue/100.);
         this->DrawShadedLimits(mode, min, max, 0.);

//       this -> DrawLegend(TString::Format("Smallest %g%",ovalue));

         break;

      // Draw a shaded band at the smallest intervals
      case 3:

         if(ovalue<50) // just to ensure there's some sense in the number
            ovalue = 68.; // default is 68%

         this -> DrawSmallest(mode,ovalue);

         break;

      // Sort out bad options and warn.
      default:

         BCLog::Out(BCLog::warning, BCLog::warning, Form("BCH1D::Draw : Invalid option %d",options));
         break;
   }
}

void BCH1D::DrawLegend ( const char *  text  ) 

Include a legend for the symbols of mean, mode, median and confidence band used in 1D marginalized posterior distributions.

Parameters:
text the text used to name the legend entry for the confidence band

Definition at line 261 of file BCH1D.cxx.

{
   //draw on top right corner

   TLegend* legend = new TLegend(0.73, 0.72, 0.86, 0.85);
   legend -> SetFillColor(kWhite);
   legend -> SetBorderSize(1);

   TMarker* triangle = new TMarker(0, 0, 23);
   triangle -> SetMarkerColor(kRed);
   legend -> AddEntry(triangle, "Global mode", "P");

   TMarker* diamond = new TMarker(0, 0, 27);
   diamond -> SetMarkerColor(kBlue);
   legend -> AddEntry(diamond, "Mean", "P");

   TLine * line;
   line = new TLine();
   line -> SetLineStyle(2);
   line -> SetLineColor(kGreen + 2);
   legend -> AddEntry(line, "Median", "l");

   TLegend* band = new TLegend(0, 0, 1, 1);
   band -> SetFillColor(kYellow);
   legend -> AddEntry(band, text, "F");

   legend -> SetTextAlign(12);

   //fine tuned by hand so text legible even with several plots in a row
   legend -> SetTextSize(0.016);

   legend -> Draw();
}

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 299 of file BCH1D.cxx.

{
   double maximum = fHistogram -> GetMaximum();

   double x0 = mode;
   double y0 = fHistogram->GetBinContent( fHistogram->FindBin(mode) );

   double x1 = fHistogram->GetMean();
   double y1 = fHistogram->GetBinContent( fHistogram->FindBin(x1) );

   double x2 = this -> GetQuantile(.5); // median
   double y2 = fHistogram->GetBinContent( fHistogram->FindBin(x2) );

   double ysize = maximum*1.3;

   double xmin = fHistogram->GetXaxis()->GetXmin();
   double xmax = fHistogram->GetXaxis()->GetXmax();

   // draw histogram with axes first
   TH2D * hsc = new TH2D(
         TString::Format("h2scale_%s_%d",fHistogram->GetName(),BCLog::GetHIndex()),"",
         10, xmin, xmax, 10, 0., ysize);
   hsc -> SetStats(0);
   hsc -> SetXTitle(fHistogram->GetXaxis()->GetTitle());
   hsc -> SetYTitle(fHistogram->GetYaxis()->GetTitle());
   hsc -> Draw();

   // draw histogram
   fHistogram -> SetLineWidth(1);
   fHistogram -> Draw("same");

   // draw yellow shaded region between min and max
   TH1D * hist_shaded = this->GetSubHisto(min,max,TString::Format("%s_sub_%d",fHistogram->GetName(),BCLog::GetHIndex()));
   hist_shaded -> SetFillStyle(1001);
   hist_shaded -> SetFillColor(kYellow);

   // draw shaded histogram
   hist_shaded -> Draw("same");

   gPad->RedrawAxis();

   // draw triangle for mode
   TPolyLine * tmax;

   double dx = 0.01*(xmax-xmin);
   double dy = 0.04*(ysize);
   y0+=dy/5.;
   double tmax_x[] = {x0, x0 + dx, x0 - dx, x0};
   double tmax_y[] = {y0, y0 + dy, y0 + dy, y0};
   tmax = new TPolyLine(4,tmax_x,tmax_y);
   tmax->SetLineColor(kRed);
   tmax->SetLineWidth(1);
   tmax->SetFillColor(kRed);
   tmax->Draw();
   tmax->Draw("f");

   // draw diamond for mean
   TPolyLine * tmean;

   y1+=dy/5.;
// double tmean_x[] = {x1, x1 + dx, x1 - dx, x1};
// double tmean_y[] = {y1, y1 + dy, y1 + dy, y1};
   double tmean_x[] = {x1, x1 + dx, x1 , x1 - dx, x1};
   double tmean_y[] = {y1, y1 + dy/2., y1 + dy, y1 + dy/2., y1};
   tmean = new TPolyLine(5,tmean_x,tmean_y);
   tmean->SetLineColor(kBlue);
// tmean->SetFillColor(kWhite);
   tmean->SetLineWidth(1);
// tmean->SetLineStyle(1);
   tmean->Draw();

/*
   // draw arrow for median
   TPolyLine * tmed;
   TPolyLine * tmed2;

// y2+=dy/5.;
   y2=0.+dy/5.;
   double tmed_x[] = {x2 + dx, x2, x2 - dx};
   double tmed_y[] = {y2 + dy, y2, y2 + dy};
   double tmed2_x[] = {x2, x2};
   double tmed2_y[] = {y2, y2 + dy*2.};
   tmed = new TPolyLine(3,tmed_x,tmed_y);
   tmed2 = new TPolyLine(2,tmed2_x,tmed2_y);
   tmed->SetLineColor(kGreen+2);
   tmed->SetLineWidth(1);
   tmed2->SetLineColor(kGreen+2);
   tmed2->SetLineWidth(1);
   tmed->Draw();
   tmed2->Draw();
*/
   // draw dashed line for median
   TLine * line;
   line = new TLine();
   line -> SetLineStyle(2);
   line -> SetLineColor(kGreen+2);
   line -> DrawLine(x2, 0., x2, y2);


   // write mode location and shaded band

   // format of the number
   double delta_max = fmax(fabs(max-x1),fabs(x1-min));

   int sd = 2 + (int)log10(fabs(x1/delta_max));

   if( (int)log10(x1) > (int)log10(delta_max) )
      sd++;

   TLatex * tmax_text = new TLatex();
   tmax_text->SetTextSize(0.045);
   tmax_text->SetTextFont(62);
   tmax_text->SetTextAlign(22); // center
// tmax_text->SetTextAlign(13); // top-left

   double xprint=(xmax+xmin)/2.;
// double xprint=xmin+(xmax-xmin)/20.;
   double yprint=ysize*(1-1.4*tmax_text->GetTextSize());

   if(fabs(limit)<50) // just to ensure there's some sense in the number
      tmax_text->DrawLatex(xprint,yprint,
         TString::Format( TString::Format("%%s^{med} = %%.%dg ^{+%%.2g}_{ -%%.2g}",sd),
            fHistogram->GetXaxis()->GetTitle(), x2, max-x2, x2-min));

   else if (limit<0)
      tmax_text->DrawLatex(xprint,yprint,
         TString::Format( TString::Format("%%s (%d%%%% prob.) < %%.4g",-(int)limit),
            fHistogram->GetXaxis()->GetTitle(), max));

   else if (limit>0)
      tmax_text->DrawLatex(xprint,yprint,
         TString::Format( TString::Format("%%s (%d%%%% prob.) > %%.4g",(int)limit),
            fHistogram->GetXaxis()->GetTitle(), min));

/*
   TLegend * leg = new TLegend(.61,.7,.9,.88);
   leg -> SetFillColor(kWhite);
   leg -> SetFillStyle(0);
   leg -> SetBorderSize(0);
   leg -> AddEntry(line,"Median", "l");
   TH1D * hh0 = new TH1D("hh0","",10,0,10);
   hh0->SetMarkerStyle(23);
   hh0->SetMarkerColor(kBlue);
   hh0->SetMarkerSize(1);
   leg -> AddEntry(hh0,"Mean","p");
   TH1D * hh1 = new TH1D("hh1","",10,0,10);
   hh1->SetMarkerStyle(23);
   hh1->SetMarkerColor(kRed);
   hh1->SetMarkerSize(1);
   leg -> AddEntry(hh1,"Mode","p");
   leg -> AddEntry(hist_shaded, "Central 68%", "f");
   leg ->Draw();
*/

}

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 457 of file BCH1D.cxx.

{
   // prepare triangle for mode
   TPolyLine * tmax;

   double x0 = mode;
   double y0 = fHistogram->GetBinContent( fHistogram->FindBin(mode) );
   double xmin = fHistogram->GetXaxis()->GetXmin();
   double xmax = fHistogram->GetXaxis()->GetXmax();
   double ysize = 1.3 * fHistogram -> GetMaximum();

   double x1 = fHistogram->GetMean();
   double y1 = fHistogram->GetBinContent( fHistogram->FindBin(x1) );

   double x2 = this -> GetQuantile(.5); // median
   double y2 = fHistogram->GetBinContent( fHistogram->FindBin(x2) );

   double dx = 0.01*(xmax-xmin);
   double dy = 0.04*(ysize);
   double tmax_x[] = {x0, x0 + dx, x0 - dx, x0};
   double tmax_y[] = {y0, y0 + dy, y0 + dy, y0};
   tmax = new TPolyLine(4,tmax_x,tmax_y);
   tmax->SetLineColor(kRed);
   tmax->SetFillColor(kRed);

   // draw histogram with axes first
   TH2D * hsc = new TH2D(
         TString::Format("h2scale_%s_%d",fHistogram->GetName(),BCLog::GetHIndex()),"",
         10, xmin, xmax, 10, 0., ysize);
   hsc -> SetStats(0);
   hsc -> SetXTitle(fHistogram->GetXaxis()->GetTitle());
   hsc -> SetYTitle(fHistogram->GetYaxis()->GetTitle());
   hsc -> Draw();

   // histogram to be filled with band
   TH1D * hist_temp1 = (TH1D*) fHistogram -> Clone();
   hist_temp1 -> Scale(1.0/fHistogram -> Integral("width"));
   hist_temp1 -> SetFillColor(kYellow);
   hist_temp1 -> SetFillStyle(1001);

   // temporary histogram
   TH1D * hist_temp2 = (TH1D*) fHistogram -> Clone();
   hist_temp2 -> Scale(1.0/fHistogram -> Integral("width"));

   // clear content
   hist_temp1 -> Reset();

   // loop over original histogram and copy bin untils integral equals
   // "prob"
   prob /= 100.;

   // integral
   double sum = 0.0;

   int n=0;
   int nbins=fHistogram->GetNbinsX();

   // loop
   while (sum < prob && n < nbins)
   {
      n++;

      // find bin with maximum
      int bin = hist_temp2 -> GetMaximumBin();

      // copy bin to new histogram
      double val = hist_temp2 -> GetBinContent(bin);
      hist_temp1 -> SetBinContent(bin, val);

      // remove maximum from temporary histogram
      hist_temp2 -> SetBinContent(bin, 0.0);

      // integrate by summing
      sum += val * hist_temp2->GetBinWidth(bin);
   }

   // scale histogram
   hist_temp1 -> Scale(fHistogram -> Integral("width"));

   // draw histograms
   fHistogram -> Draw("same");
   hist_temp1 -> Draw("same");

   // draw triangle for mode
   tmax->Draw("f");

   if(drawmean)
   {
      // draw triangle for mean
      // draw diamond for mean
      TPolyLine * tmean;

      y1+=dy/5.;
//    double tmean_x[] = {x1, x1 + dx, x1 - dx, x1};
//    double tmean_y[] = {y1, y1 + dy, y1 + dy, y1};
      double tmean_x[] = {x1, x1 + dx, x1 , x1 - dx, x1};
      double tmean_y[] = {y1, y1 + dy/2., y1 + dy, y1 + dy/2., y1};
      tmean = new TPolyLine(5,tmean_x,tmean_y);
      tmean->SetLineColor(kBlue);
//    tmean->SetFillColor(kWhite);
      tmean->SetLineWidth(1);
//    tmean->SetLineStyle(1);
      tmean->Draw();

      // draw dashed line for median
      TLine * line;
      line = new TLine();
      line -> SetLineStyle(2);
      line -> SetLineColor(kGreen+2);
      line -> DrawLine(x2, 0., x2, y2);
   }

   // free memory
   delete hist_temp2;
}

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

Definition at line 57 of file BCH1D.h.

         { 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.

{
   double integral = 0;

   int binmin = fHistogram -> FindBin(valuemin);
   int binmax = fHistogram -> FindBin(valuemax);

   // use ROOT function to calculate integral.
   integral = fHistogram -> Integral(binmin, binmax);

   return integral;
}

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.

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

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

Definition at line 62 of file BCH1D.h.

         { return fHistogram -> GetMean(); };

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

Definition at line 71 of file BCH1D.h.

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

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

Definition at line 48 of file BCH1D.cxx.

{
   return fHistogram -> GetBinCenter(fHistogram -> GetMaximumBin());
}

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.

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

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.

{
   int nquantiles = 1;
   double quantiles[1];
   double probsum[1];

   probsum[0] = probability;

   // use ROOT function to calculat quantile.
   fHistogram -> GetQuantiles(nquantiles, quantiles, probsum);

   return quantiles[0];
}

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

Definition at line 91 of file BCH1D.h.

         { 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 575 of file BCH1D.cxx.

{
   if(content<=0. || content>= 1.)
      return 0.;

   int nbins = fHistogram -> GetNbinsX();

   double factor = fHistogram->Integral("width");
   if(factor==0)
      return 0.;

   // normalize if not yet done
   fHistogram->Scale(1./factor);

   double xmin = fHistogram->GetXaxis()->GetXmin();
   double xmax = fHistogram->GetXaxis()->GetXmax();
   double width = xmax - xmin;

   double xdown=xmin;
   double xup=xmax;

   int ndiv = 100;
   if(nbins<100)
      ndiv = 1000;
   if(nbins>1000)
      ndiv = 10;

   int warn=0;

   double integral_out=0.;

   // loop through the bins
   for(int i=1;i<nbins+1;i++)
   {
      if(fHistogram->Integral(i,nbins,"width") < content)
         break;

      // get width of the starting bin
      double firstbinwidth = fHistogram->GetBinWidth(i);

      // divide i-th bin in ndiv divisions
      // integrate starting at each of these divisions
      for(int j=0;j<ndiv;j++)
      {
         double dxdown = (double)(ndiv-j)/(double)ndiv * firstbinwidth;
         xdown = fHistogram->GetBinLowEdge(i) + firstbinwidth - dxdown;
         double integral = dxdown * fHistogram->GetBinContent(i);

         if(integral>content)
         {
            // content fits within one bin
            xup = xdown + content / fHistogram->GetBinContent(i);
            warn = 1;
         }
         else
            for(int k=i+1;k<nbins+1;k++)
            {

               double thisbin = fHistogram->GetBinContent(k) * fHistogram->GetBinWidth(k);

               if(integral+thisbin==content)
               {
                  xup = fHistogram->GetBinLowEdge(k+1);
                  break;
               }

               if(integral+thisbin>content)
               {
                  xup = fHistogram->GetBinLowEdge(k) + (content-integral)/thisbin * fHistogram->GetBinWidth(k);
                  integral += thisbin * (xup-fHistogram->GetBinLowEdge(k))/fHistogram->GetBinWidth(k);
                  break;
               }

               integral += thisbin;
            }

         if(integral < content)
            continue;

         if(xup - xdown < width)
         {
            // store necessary information
            width = xup - xdown;
            xmin  = xdown;
            xmax  = xup;
            integral_out = integral;
         }
      }
   }

   if(warn)
   {
      BCLog::Out(BCLog::warning,BCLog::warning,
            Form("BCH1D::GetSmallestInterval : The requested content of %d%% fits within one bin.",(int)(content*100)));
      BCLog::Out(BCLog::warning,BCLog::warning,
            "BCH1D::GetSmallestInterval : MAKE FINER BINNING (OR CHANGE RANGE) !!!");
   }

   // restore normalization to state before calling this method
   fHistogram->Scale(factor);

   min=xmin;
   max=xmax;

   return integral_out;
}

TH1D * BCH1D::GetSmallestIntervalHistogram ( double  level  ) 

Definition at line 802 of file BCH1D.cxx.

{
   // create a new histogram which will be returned and all yellow
   TH1D * hist_yellow = (TH1D*) fHistogram -> Clone();
   hist_yellow -> Reset();
   hist_yellow -> SetFillStyle(1001);
   hist_yellow -> SetFillColor(kYellow);

   // copy a temporary histogram
   TH1D * hist_temp = (TH1D*) fHistogram -> Clone(TString::Format("%s_%i",fHistogram->GetName(),BCLog::GetHIndex()));
   double factor = hist_temp -> Integral("");

   if(factor == 0)
      return 0;

   hist_temp -> Scale(1. / factor);

   // here's the algorithm:
   // 1. find the maximum bin in the temporary histogram and copy it to
   // the yellow histogram.
   // 2. remove this bin from the temporary histogram.
   // 3. repeat this until a total of "level" probability is copied to
   // the yellow histogram.

   double sumprob = 0.0;

   while (sumprob < level)
   {
      // find maximum bin and value
      int bin = hist_temp -> GetMaximumBin();
      double value = hist_temp -> GetMaximum();

      // copy "1" into the corresponding bin in the yellow histogram
      hist_yellow -> SetBinContent(bin, 1.0);

      // set the bin value in the temporary histogram to zero
      hist_temp -> SetBinContent(bin, 0.0);

      // increase probability sum
      sumprob += value;
   }

   delete hist_temp;

   return hist_yellow;
}

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

Definition at line 851 of file BCH1D.cxx.

{
   std::vector <double> v;

   TH1D * hist = this -> GetSmallestIntervalHistogram(content);

   int nbins = hist -> GetNbinsX();
   int ninter = 0;
   int lastbin = -1;

   double max = -1;
   double localmax = -1;
   double localmaxpos = -1;
   double localint = 0;
   bool flag = false;

   for (int i = 1; i <= nbins; ++i)
   {
      // interval starts here
      if (!flag && hist -> GetBinContent(i) > 0.)
      {
         flag = true;
         v.push_back(hist -> GetBinLowEdge(i));

         // remember start position of the interval
         lastbin = i;

         // increase number of intervals
         ninter++;

         // reset local maximum
         localmax = fHistogram -> GetBinContent(i);
         localmaxpos = hist -> GetBinLowEdge(i);

         // reset local integral
         localint = 0;
      }

      // interval stops here
      if ((flag && !(hist -> GetBinContent(i) > 0.)) || (flag && i == nbins))
      {
         flag = false;
         v.push_back(hist -> GetBinLowEdge(i) + hist -> GetBinWidth(i));

         // set right bin to maximum if on edge
         if (i == nbins && localmax < fHistogram -> GetBinContent(i))
            localmaxpos = hist -> GetBinCenter(i) + 0.5 * hist -> GetBinWidth(i);

         // find the absolute maximum
         if (localmax > max)
            max = localmax;

         // save local maximum
         v.push_back(localmax);
         v.push_back(localmaxpos);

         // save local integral
         v.push_back(localint);
      }

      // find local maximum
      if (i < nbins && localmax < fHistogram -> GetBinContent(i))
      {
         localmax = fHistogram -> GetBinContent(i);
         localmaxpos = hist -> GetBinCenter(i);
      }

      // increase area
      localint += fHistogram -> GetBinContent(i) / fHistogram -> Integral();
   }

   // rescale absolute heights to relative heights
   for (int i = 0; i < ninter; ++i)
      v[i*5+2] = v.at(i*5+2) / max;

   return v;
}

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 723 of file BCH1D.cxx.

{
   if(min>max)
   {
      double t=min;
      min=max;
      max=t;
   }

   int imin = fHistogram->FindBin(min);
   int imax = fHistogram->FindBin(max);

   int nbins = fHistogram->GetNbinsX();
   double xmin = fHistogram->GetXaxis()->GetXmin();
   double xmax = fHistogram->GetXaxis()->GetXmax();

   if( min==max || (min<=xmin && max>=xmax) )
   {
      TH1D * h0 = (TH1D*) fHistogram->Clone(name);
      return h0;
   }

   if (imin<1)
   {
      imin=1;
      min=xmin;
   }
   if (imax>nbins)
   {
      imax=nbins;
      max=xmax;
   }

   double * xb = new double[nbins+3]; // nbins+1 original bin edges + 2 new bins
   int n=0; // counter

   int domin=1;
   int domax=1;

   for (int i=1;i<nbins+2;i++)
   {
      double x0 = fHistogram->GetBinLowEdge(i);

      if (min<x0 && domin)
      {
         xb[n++]=min;
         domin=0;
      }
      else if (min==x0)
         domin=0;

      if (max<x0 && domax)
      {
         xb[n++]=max;
         domax=0;
      }
      else if (max==x0)
         domax=0;

      xb[n++]=x0;
   }

   // now define the new histogram
   TH1D * h0 = new TH1D(name,"",n-1,xb);
   for(int i=1;i<n;i++)
   {
      double x0 = h0->GetBinCenter(i);
      if(x0<min || x0>max)
         continue;

      int bin=fHistogram->FindBin(x0);
      h0->SetBinContent(i, fHistogram->GetBinContent(bin));
   }

   return h0;
}

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 684 of file BCH1D.cxx.

{
   int imin = fHistogram->FindBin(min);
   int imax = fHistogram->FindBin(max);

   int nbins = fHistogram->GetNbinsX();

   // if outside of histogram range, return -1.
   if ( imin<1 || imin>nbins || imax<1 || imax>nbins )
      return -1.;

   if ( imin==imax )
      return -1.;

   // swap values if necessary
   if (imin>imax)
   {
      int i=imin;
      double x=min;
      imin=imax, min=max;
      imax=i, max=x;
   }

   // calculate first bin
   double first = ( fHistogram->GetBinLowEdge(imin+1) - min ) * fHistogram->GetBinContent(imin);

   // calculate last bin
   double last = ( max - fHistogram->GetBinLowEdge(imax) ) * fHistogram->GetBinContent(imax);

   // calculate inbetween
   double inbetween=0.;
   if(imax-imin>1)
      inbetween = fHistogram->Integral(imin+1, imax-1, "width");

   return first + inbetween + last;
}

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.

{
   char file[256];
   int i=0;
   while(i<255 && *filename!='\0')
      file[i++]=*filename++;
   file[i]='\0';

// temporary
   fHistogram -> SetLineWidth(1);

   // create temporary canvas
   TCanvas * c;
   if(ww > 0 && wh > 0)
      c = new TCanvas("c","c",ww,wh);
   else
      c = new TCanvas("c","c");

   c -> cd();
   this->Draw(options, ovalue);

   gPad->RedrawAxis();

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

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.

{
   // check if limit is between 68% and 100%. Give a warning if not ...
   if(limit>=100. || limit<68.)
      BCLog::Out(BCLog::warning,BCLog::warning,
         Form("BCH1D::SetDefaultCLLimit : Value %.1f out of range. Keeping default %.1f%% CL limit.",limit,fDefaultCLLimit));

   // ... or set value if true.
   else
      fDefaultCLLimit = limit;
}

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

Set global mode

Definition at line 125 of file BCH1D.h.

         { fMode=mode; fModeFlag=1; };

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

Sets the histogram.

Definition at line 115 of file BCH1D.h.

         { fHistogram = hist; };


Member Data Documentation

double BCH1D::fDefaultCLLimit [private]

Default confidence level limit

Definition at line 221 of file BCH1D.h.

TH1D* BCH1D::fHistogram [private]

The 1D histogram

Definition at line 217 of file BCH1D.h.

double BCH1D::fMode [private]

Global mode

Definition at line 225 of file BCH1D.h.

int BCH1D::fModeFlag [private]

"Is there a global mode?" flag

Definition at line 229 of file BCH1D.h.


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