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

#include <BCH1D.h>

Public Member Functions
Constructors and destructors
	BCH1D (TH1D *hist=0)

	~BCH1D ()

Member functions (get)
TH1D *	GetHistogram ()

double	GetMean ()

double	GetMode ()

double	GetMedian ()

double	GetQuantile (double probablity)

double	GetLimit (double probability)

double	GetRMS ()

double	GetSTD ()

double	GetVariance ()

double	GetSkew ()

double	GetKurtosis ()

double	GetIntegral (double valuemin, double valuemax)

double	GetPValue (double probability)

int	GetColor (int index)

Member functions (set)
void	SetColorScheme (int scheme)

void	SetHistogram (TH1D *hist)

void	SetDefaultCLLimit (double limit)

void	SetGlobalMode (double mode)

Member functions (miscellaneous methods)
void	Print (const char *filename, std::string options="BTsiB3CS1D0Lmeanmode", std::vector< double > intervals=std::vector< double >(0), int ww=0, int wh=0)

void	Print (const char *filename, std::string options, double interval, int ww=0, int wh=0)

void	Draw (std::string options="BTsiB3CS1D0Lmeanmode", std::vector< double > intervals=std::vector< double >(0))

void	Draw (std::string options, double interval)

void	DrawDelta (double value) const

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)

Static Private Member Functions
static unsigned int	getNextIndex ()

Private Attributes
TH1D *	fHistogram

double	fDefaultCLLimit

double	fMode

int	fModeFlag

std::vector< int >	fColors

std::vector< TObject * >	fROOTObjects

Static Private Attributes
static unsigned int	fHCounter =0

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 31 of file BCH1D.h.

Constructor & Destructor Documentation

BCH1D::BCH1D ( TH1D * hist = 0)

The default constructor.

Definition at line 37 of file BCH1D.cxx.

   : fHistogram(hist)
   , fDefaultCLLimit(95.) // in percent
   , fMode(0)
   , fModeFlag(0)
 {}

BCH1D::~BCH1D ( )

The default destructor.

Definition at line 45 of file BCH1D.cxx.

 {
    for (unsigned i = 0; i < fROOTObjects.size(); ++i)
       delete fROOTObjects[i];
 }

Member Function Documentation

void BCH1D::Draw	(	std::string	options = `"BTsiB3CS1D0Lmeanmode"`,
		std::vector< double >	intervals = `std::vector<double>(0)`
	)

Draw distribution into the active canvas.

Parameters

options Drawing options:
BTci : band type is central interval [default]
BTsi : band type is/are smallest interval(s)
BTul : band type is upper limit
BTll : band type is lower limit
B1 : draw one band between values specified in intervals [default]
B2 : draw two bands between values specified in intervals
B3 : draw three bands between values specified in intervals
D0 : draw histogram [default]
D1 : draw smooth curve
CS0 : choose color scheme 0 (B&W)
CS1 : choose color scheme 1 (green/yellow/red) [default]
CS2 : choose color scheme 2 (blueish colors)
CS3 : choose color scheme 3 (redish colors)
smooth1 : use ROOT smoothing algorithm once
smooth3 : use ROOT smoothing algorithm three times
smooth5 : use ROOT smoothing algorithm five times
smooth10 : use ROOT smoothing algorithm ten times
median : draw median and central interval
mode : draw global mode and standard deviation
mean : draw mean value and standard deviation
quartiles : indicate quartiles
deciles : indicate deciles
percentiles : indicate percentiles
L : add a legend
same: add histogram on top of another histogram

intervals the intervals

Definition at line 231 of file BCH1D.cxx.

 {
    // option flags
    bool flag_legend = false;
    bool flag_logy = false;
    bool flag_mode = false;
    bool flag_median = false;
    bool flag_mean = false;
    bool flag_quartiles = false;
    bool flag_deciles = false;
    bool flag_percentiles = false;
    bool flag_smooth1 = false;
    bool flag_smooth3 = false;
    bool flag_smooth5 = false;
    bool flag_smooth10 = false;
 
    // band type
    int bandtype = 0;
 
    // number of bands
    int nbands = 0; // number of shaded bands
 
    // define draw options called in TH1D::Draw(...)
    std::string draw_options = "";
 
    // check content of options string
    if (options.find("smooth1") < options.size()) {
       flag_smooth1 = true;
    }
 
    if (options.find("smooth3") < options.size()) {
       flag_smooth3 = true;
    }
 
    if (options.find("smooth5") < options.size()) {
       flag_smooth5 = true;
    }
 
    if (options.find("smooth10") < options.size()) {
       flag_smooth10 = true;
    }
 
    if (options.find("same") < options.size()) {
       draw_options.append("SAME");
    }
 
    if (options.find("logy") < options.size()) {
       flag_logy = true;
    }
 
    if (options.find("L") < options.size()) {
       flag_legend = true;
    }
 
    if (options.find("D1") < options.size()) {
       draw_options.append("C");
    }
 
    if (options.find("BTsi") < options.size()) {
       bandtype = 1;
    }
    else if (options.find("BTul") < options.size()) {
       bandtype = 2;
    }
    else if (options.find("BTll") < options.size()) {
       bandtype = 3;
    }
    else {
       bandtype = 0;
    }
 
    if (options.find("B1") < options.size()) {
       nbands = 1;
       if (bandtype == 0 && intervals.size() != 2) {
          intervals.clear();
          intervals.push_back(0.1587);
          intervals.push_back(0.8413);
       }
       else if (bandtype == 1 && intervals.size() != 1) {
          intervals.clear();
          intervals.push_back(0.6827);
       }
       else if (bandtype == 2 && intervals.size() != 1) {
          intervals.clear();
          intervals.push_back(0.90);
       }
       else if (bandtype == 3 && intervals.size() != 1) {
          intervals.clear();
          intervals.push_back(0.10);
       }
    }
 
    if (options.find("B2") < options.size()) {
       nbands = 2;
       if (bandtype == 0 && intervals.size() != 4) {
          intervals.clear();
          intervals.push_back(0.1587);
          intervals.push_back(0.8413);
          intervals.push_back(0.0228);
          intervals.push_back(0.9772);
       }
       else if (bandtype == 1 && intervals.size() != 2) {
          intervals.clear();
          intervals.push_back(0.6827);
          intervals.push_back(0.9545);
       }
       else if (bandtype == 2 && intervals.size() != 2) {
          intervals.clear();
          intervals.push_back(0.90);
          intervals.push_back(0.95);
       }
       else if (bandtype == 3 && intervals.size() != 2) {
          intervals.clear();
          intervals.push_back(0.10);
          intervals.push_back(0.05);
       }
    }
 
    if (options.find("B3") < options.size()) {
       nbands = 3;
       if (bandtype == 0 && intervals.size() != 6) {
          intervals.clear();
          intervals.push_back(0.1587);
          intervals.push_back(0.8413);
          intervals.push_back(0.0228);
          intervals.push_back(0.9772);
          intervals.push_back(0.0013);
          intervals.push_back(0.9987);
       }
       else if (bandtype == 1 && intervals.size() != 3) {
          intervals.clear();
          intervals.push_back(0.6827);
          intervals.push_back(0.9545);
          intervals.push_back(0.9973);
       }
       else if (bandtype == 2 && intervals.size() != 3) {
          intervals.clear();
          intervals.push_back(0.90);
          intervals.push_back(0.95);
          intervals.push_back(0.99);
       }
       else if (bandtype == 3 && intervals.size() != 3) {
          intervals.clear();
          intervals.push_back(0.10);
          intervals.push_back(0.05);
          intervals.push_back(0.01);
       }
    }
 
    if (options.find("CS0") < options.size()) {
       SetColorScheme(0);
    }
    else if (options.find("CS1") < options.size()) {
       SetColorScheme(1);
    }
    else if (options.find("CS2") < options.size()) {
       SetColorScheme(2);
    }
    else if (options.find("CS3") < options.size()) {
       SetColorScheme(3);
    }
    else  {
       SetColorScheme(1);
    }
 
    if (options.find("mode") < options.size()) {
       if (fModeFlag)
          flag_mode = true;
    }
    if (options.find("median") < options.size()) {
       flag_median = true;
    }
    if (options.find("mean") < options.size()) {
       flag_mean = true;
    }
    if (options.find("quartiles") < options.size()) {
       flag_quartiles = true;
    }
    if (options.find("deciles") < options.size()) {
       flag_deciles = true;
    }
    if (options.find("percentiles") < options.size()) {
       flag_percentiles = true;
    }
 
    // normalize histogram to unity
    fHistogram->Scale(1./fHistogram->Integral("width"));
 
    // prepare legend
    TLegend* legend = new TLegend();
    legend->SetBorderSize(0);
    legend->SetFillColor(kWhite);
    legend->SetTextAlign(12);
    legend->SetTextFont(62);
    legend->SetTextSize(0.03);
 
    // add legend to list of objects
    fROOTObjects.push_back(legend);
 
    // smooth
    if (flag_smooth1) {
       fHistogram->Smooth(1);
       fHistogram->Scale(1.0/fHistogram->Integral("width"));
    }
    if (flag_smooth3) {
       fHistogram->Smooth(3);
       fHistogram->Scale(1.0/fHistogram->Integral("width"));
    }
 
    if (flag_smooth5) {
       fHistogram->Smooth(5);
       fHistogram->Scale(1.0/fHistogram->Integral("width"));
    }
    if (flag_smooth10) {
       fHistogram->Smooth(10);
       fHistogram->Scale(1.0/fHistogram->Integral("width"));
    }
 
    // draw histogram
    fHistogram->Draw(draw_options.c_str());
 
    // draw bands
    for (int i = 0; i < nbands; ++i) {
       int col = GetColor(nbands-i-1);
 
       double prob_low  = 0;
       double prob_high = 0;
       double prob_interval = 0;
       double xlow  = 0;
       double xhigh = 0;
 
       TH1D * hist_band = 0;
 
       if (bandtype == 0) {
          prob_low  = intervals[2*(nbands-i)-2];
          prob_high = intervals[2*(nbands-i)-1];
          xlow  = GetQuantile(prob_low);
          xhigh = GetQuantile(prob_high);
          prob_interval = prob_high - prob_low;
 
          hist_band = GetSubHisto(xlow, xhigh, TString::Format("%s_sub_%d", fHistogram->GetName(), BCLog::GetHIndex()));
       }
       else if (bandtype == 1) {
          prob_interval = GetSmallestInterval(xlow, xhigh, intervals[nbands-1-i]);
          hist_band = GetSmallestIntervalHistogram(intervals[nbands-1-i]);
          for (int ibin = 1; ibin < hist_band->GetNbinsX(); ++ibin)
             hist_band->SetBinContent(ibin, hist_band->GetBinContent(ibin)*fHistogram->GetBinContent(ibin));
       }
       else if(bandtype == 2) {
          xlow = 0.;
          xhigh = GetQuantile(intervals[nbands-1-i]);
          hist_band = GetSubHisto(xlow, xhigh, TString::Format("%s_sub_%d", fHistogram->GetName(), BCLog::GetHIndex()));
          prob_interval = intervals[nbands-1-i];
       }
       else if(bandtype == 3) {
          xlow = GetQuantile(intervals[nbands-1-i]);
          xhigh = GetQuantile(1.);
          hist_band = GetSubHisto(xlow, xhigh, TString::Format("%s_sub_%d", fHistogram->GetName(), BCLog::GetHIndex()));
          prob_interval = 1.-intervals[nbands-1-i];
       }
 
       // set style of band histogram
       hist_band->SetFillStyle(1001);
       hist_band->SetFillColor(col);
       hist_band->SetLineColor(col);
 
       // draw shaded histogram
       hist_band->Draw(std::string(draw_options+std::string("same")).c_str());
 
       // draw histogram again
       fHistogram->Draw(std::string(std::string("SAME")+draw_options).c_str());
 
       // add to legend
       std::string legend_label;
       if (bandtype == 0)
          legend_label.append(Form("central %.1f%% interval ", prob_interval*100));
       else if (bandtype == 1)
          legend_label.append(Form("smallest %.1f%% interval(s)", prob_interval*100));
       else if (bandtype == 2)
          legend_label.append(Form("%.0f%% upper limit", prob_interval*100));
       else if (bandtype == 3)
          legend_label.append(Form("%.0f%% lower limit", prob_interval*100));
 
       legend->AddEntry(hist_band, legend_label.c_str(), "F");
 
       // add hist_band to list of objects
       fROOTObjects.push_back(hist_band);
    }
 
    gPad->RedrawAxis();
 
    // prepare size of histogram
    double ymin     = 0;
    double ymaxhist = fHistogram->GetBinContent(fHistogram->GetMaximumBin());
    double ymax     = ymaxhist;
 
    // check if log axis
    if (flag_logy)
       ymin = 1e-4*ymaxhist;
 
    // quantiles
    TLine* line_quantiles = new TLine();
    line_quantiles->SetLineStyle(2);
    line_quantiles->SetLineColor(GetColor(3));
 
    if (flag_quartiles) {
       for (int i = 1; i < 4; ++i) {
          double quantile_x = GetQuantile(0.25*i);
          int quantile_xbin = fHistogram->FindBin(quantile_x);
          double quantile_y = fHistogram->GetBinContent(quantile_xbin);
          double quantile_ymin = 0;
          if (flag_logy)
             quantile_ymin = 1e-4*ymaxhist;
          line_quantiles->DrawLine(quantile_x, quantile_ymin, quantile_x, quantile_y);
       }
       TLegendEntry* le = legend->AddEntry(line_quantiles, "quartiles", "L");
       if (nbands>0)
          le->SetFillColor(GetColor(0));
       else
          le->SetFillColor(kWhite);
       le->SetFillStyle(1001);
    }
 
    if (flag_deciles) {
       for (int i = 1; i < 10; ++i) {
          double quantile_x = GetQuantile(0.10*i);
          int quantile_xbin = fHistogram->FindBin(quantile_x);
          double quantile_y = fHistogram->GetBinContent(quantile_xbin);
          double quantile_ymin = 0;
          if (flag_logy)
             quantile_ymin = 1e-4*ymaxhist;
          line_quantiles->DrawLine(quantile_x, quantile_ymin, quantile_x, quantile_y);
       }
       TLegendEntry* le = legend->AddEntry(line_quantiles, "deciles", "FL");
       le->SetFillColor(GetColor(0));
       le->SetFillStyle(1001);
    }
 
    if (flag_percentiles) {
       for (int i = 1; i < 100; ++i) {
          double quantile_x = GetQuantile(0.01*i);
          int quantile_xbin = fHistogram->FindBin(quantile_x);
          double quantile_y = fHistogram->GetBinContent(quantile_xbin);
          double quantile_ymin = 0;
          if (flag_logy)
             quantile_ymin = 1e-4*ymaxhist;
          line_quantiles->DrawLine(quantile_x, quantile_ymin, quantile_x, quantile_y);
       }
       TLegendEntry* le = legend->AddEntry(line_quantiles, "percentiles", "L");
       le->SetFillColor(GetColor(0));
       le->SetFillStyle(1001);
    }
 
   // add line_quantiles to list of ROOT objects
   fROOTObjects.push_back(line_quantiles);
 
   // mean, mode, median
   TMarker* marker_mode = new TMarker(fMode, 0.50*ymaxhist, 24);
   marker_mode->SetMarkerColor(GetColor(4));
   marker_mode->SetMarkerSize(1.5);
 
   TMarker* marker_mean = new TMarker(GetMean(), 0.55*ymaxhist, 20);
   marker_mean->SetMarkerColor(GetColor(4));
   marker_mean->SetMarkerSize(1.5);
 
   TMarker* marker_median = new TMarker(GetMedian(), 0.45*ymaxhist, 21);
   marker_median->SetMarkerColor(GetColor(4));
   marker_median->SetMarkerSize(1.5);
 
   // mode
   TArrow* arrow_mode = new TArrow(fMode, 0.485*ymaxhist,
               fMode, ymin,
               0.02, "|>");
   arrow_mode->SetLineColor(GetColor(4));
   arrow_mode->SetFillColor(GetColor(4));
 
   // standard deviation
   TArrow* arrow_std = new TArrow(GetMean()-GetRMS(), 0.55*ymaxhist,
              GetMean()+GetRMS(), 0.55*ymaxhist,
              0.02, "<|>");
   arrow_std->SetLineColor(GetColor(4));
   arrow_std->SetFillColor(GetColor(4));
 
   // central interval
   TArrow* arrow_ci = new TArrow(GetQuantile(0.1587), 0.45*ymaxhist,
             GetQuantile(0.8413), 0.45*ymaxhist,
             0.02, "<|>");
   arrow_ci->SetLineColor(GetColor(4));
   arrow_ci->SetFillColor(GetColor(4));
 
   // add marker_mean and arrow_std to list of ROOT objects
   fROOTObjects.push_back(marker_mode);
   fROOTObjects.push_back(marker_mean);
   fROOTObjects.push_back(marker_median);
   fROOTObjects.push_back(arrow_std);
   fROOTObjects.push_back(arrow_ci);
 
   if (flag_mode) {
     arrow_mode->Draw();
     marker_mode->Draw();
     TLegendEntry* le = legend->AddEntry(marker_mode, "global mode", "P");
     le->SetMarkerStyle(24);
     le->SetMarkerSize(1.5);
     le->SetMarkerColor(GetColor(4));
   }
 
   if (flag_mean) {
     arrow_std->Draw();
     marker_mean->Draw();
     TLegendEntry* le = legend->AddEntry(arrow_std, "mean and standard deviation", "PL");
     le->SetLineColor(GetColor(4));
     le->SetMarkerStyle(20);
     le->SetMarkerSize(1.5);
     le->SetMarkerColor(GetColor(4));
   }
 
   if (flag_median) {
     arrow_ci->Draw();
     marker_median->Draw();
     TLegendEntry* le = legend->AddEntry(arrow_ci, "median and central 68.3% interval", "PL");
     le->SetLineColor(GetColor(4));
     le->SetMarkerStyle(21);
     le->SetMarkerSize(1.5);
     le->SetMarkerColor(GetColor(4));
   }
 
   // calculate legend height in NDC coordinates
   double height = 0.03*legend->GetNRows();
 
   // make room for legend
   if (flag_legend)
     ymax*=(1.15+height);
   else
     ymax*=1.1;
 
   fHistogram->GetYaxis()->SetRangeUser(ymin, 1.05*ymaxhist);
 
   // calculate dimensions in NDC variables
 
   if (flag_legend)
     gPad->SetTopMargin(0.02);
 
    double xlegend1 = gPad->GetLeftMargin() + 0.10 * (1.0 - gPad->GetRightMargin() - gPad->GetLeftMargin());
    double xlegend2 = 1.0-gPad->GetRightMargin();
    double ylegend1 = 1.-gPad->GetTopMargin()-height;
    double ylegend2 = 1.-gPad->GetTopMargin();
 
    // place legend on top of histogram
    legend->SetX1NDC(xlegend1);
    legend->SetX2NDC(xlegend2);
    legend->SetY1NDC(ylegend1);
    legend->SetY2NDC(ylegend2);
 
    // draw legend
    if (flag_legend) {
       legend->Draw();
    }
 
    // rescale top margin
    gPad->SetTopMargin(1.-ylegend1+0.01);
 
    gPad->RedrawAxis();
 
    return;
 }

void BCH1D::Draw	(	std::string	options,
		double	interval
	)

Draw distribution into the active canvas.

Parameters

options Drawing options,

See Also: Print(const char * filename, std::string options, double interval, int ww=0, int wh=0)

Parameters

interval an upper or lower limit

Definition at line 698 of file BCH1D.cxx.

 {
    std::vector<double> tempvec;
    tempvec.push_back(interval);
    Draw(options, tempvec);
 }

void BCH1D::DrawDelta ( double value) const

Draw the 1D marginal for a parameter fixed by a delta prior.

Parameters

value The fixed value of the parameter.

Definition at line 706 of file BCH1D.cxx.

 {
    // draw histogram with axes first
    double xmin = fHistogram->GetXaxis()->GetXmin();
    double xmax = fHistogram->GetXaxis()->GetXmax();
    double ysize = 1.3 * fHistogram->GetMaximum();
 
    TH2D*  hsc = new TH2D(
          TString::Format("h2scale_%s_%d",fHistogram->GetName(),BCLog::GetHIndex()),"",
          50, xmin, xmax, 10, 0., ysize);
    hsc->SetStats(0);
    hsc->SetXTitle(fHistogram->GetXaxis()->GetTitle());
    hsc->SetYTitle(fHistogram->GetYaxis()->GetTitle());
    hsc->Draw();
 
    // write mode location
 
    TLatex * tmax_text = new TLatex();
    tmax_text->SetTextSize(0.035);
    tmax_text->SetTextFont(62);
    tmax_text->SetTextAlign(22); // center
 
    double xprint=(xmax+xmin)/2.;
    double yprint=ysize*(1-1.4*tmax_text->GetTextSize());
 
    tmax_text->DrawLatex(xprint,yprint, TString::Format("%s = %g", fHistogram->GetXaxis()->GetTitle(), value));
    delete tmax_text;
 
    // create a temporary histogram, to hold only one non-zero bin
    TH1D * hist_temp = new TH1D(TString::Format("h1scale_%s_%d", fHistogram->GetName(), BCLog::GetHIndex()), "", 100, xmin, xmax);
    hist_temp->SetBinContent(hist_temp->FindBin(value), fHistogram->GetMaximum());
    hist_temp->Draw("same");
    fROOTObjects.push_back(hist_temp);
 }

int BCH1D::GetColor ( int index)

inline

Returns a color of the current color scheme.

Parameters

index the color index

Returns: the color number.

Definition at line 128 of file BCH1D.h.

129 { return fColors.at(index); };

TH1D* BCH1D::GetHistogram ( )

inline

Returns: The one-dimensional histogram.

Definition at line 53 of file BCH1D.h.

54 { 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 73 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::GetKurtosis ( )

inline

Returns: The STD of the distribution.

Definition at line 107 of file BCH1D.h.

108 { return fHistogram->GetKurtosis(); };

double BCH1D::GetLimit ( double probability)

inline

Return the quantile of the distribution

Parameters

probability The probability.

Returns: The quantile of the distribution for the probability.

See Also: GetQuantile(double probablity)

Definition at line 82 of file BCH1D.h.

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

double BCH1D::GetMean ( )

inline

Returns: The mean of the distribution.

Definition at line 58 of file BCH1D.h.

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

double BCH1D::GetMedian ( )

inline

Returns: The median of the distribution.

Definition at line 67 of file BCH1D.h.

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

double BCH1D::GetMode ( )

Returns: The mode of the distribution.

Definition at line 52 of file BCH1D.cxx.

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

static unsigned int BCH1D::getNextIndex ( )

inlinestaticprivate

Helper method to get an unique number to be used in histogram naming

Definition at line 315 of file BCH1D.h.

316 { return ++fHCounter; }

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

 {
    int nquantiles = 1;
    double quantiles[1];
    double probsum[1];
 
    probsum[0] = probability;
 
    // use ROOT function to calculate quantile.
    fHistogram->GetQuantiles(nquantiles, quantiles, probsum);
 
    return quantiles[0];
 }

double BCH1D::GetRMS ( )

inline

Returns: The RMS of the distribution.

Definition at line 87 of file BCH1D.h.

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

double BCH1D::GetSkew ( )

inline

Returns: The skew of the distribution.

Definition at line 102 of file BCH1D.h.

103 { return fHistogram->GetSkewness(); };

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 742 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 = 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::OutWarning(
             Form("BCH1D::GetSmallestInterval : The requested content of %d%% fits within one bin.",(int)(content*100)));
       BCLog::OutWarning("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)

Create a histogram with the smallest intervals of the original histogram containing a certain level. The histogram is yellow.

Parameters

level the level or content of the histogram

Returns: the histogram.

Definition at line 960 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)

Return a vector containing information about the set of smallest intervals.
0 : x_min
1 : x_max
2 : relative height 3 : local mode
4 : relative area.

Parameters

content The content of the smallest interval

Returns: the vector.

Definition at line 1008 of file BCH1D.cxx.

 {
    std::vector<double> v;
 
    TH1D * hist = GetSmallestIntervalHistogram(content);
 
    int nbins = hist->GetNbinsX();
    int ninter = 0;
 
    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));
 
          // 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;
 
    delete hist;
 
    return v;
 }

double BCH1D::GetSTD ( )

inline

Returns: The standard deviation of the distribution.

Definition at line 92 of file BCH1D.h.

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

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 882 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;
    }
 
    std::vector<double> xb(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[0]);
    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::GetVariance ( )

inline

Returns: The variance of the distribution.

Definition at line 97 of file BCH1D.h.

98 { return (GetSTD()*GetSTD()); };

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 844 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,
		std::string	options = `"BTsiB3CS1D0Lmeanmode"`,
		std::vector< double >	intervals = `std::vector<double>(0)`,
		int	ww = `0`,
		int	wh = `0`
	)

Print distribution into a PostScript file.

Parameters

filename	Output filename
option	the draw options (see Draw()), plus logx : draw x-axis in log-scale logy : draw y-axis in log-scale R : rescale canvas to have a squared histogram
intervals	the intervals for the bands
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 150 of file BCH1D.cxx.

 {
    char file[256];
    int i=0;
    while(i<255 && *filename!='\0')
       file[i++]=*filename++;
    file[i]='\0';
 
    // option flags
    bool flag_logx = false;
    bool flag_logy = false;
    bool flag_rescale = false;
 
    // check content of options string
    if (options.find("logx") < options.size()) {
       flag_logx = true;
    }
 
    if (options.find("logy") < options.size()) {
       flag_logy = true;
    }
 
    if (options.find("R") < options.size()) {
       flag_rescale = true;
    }
 
    // create temporary canvas
    TCanvas * c;
    unsigned int cindex = getNextIndex();
    if(ww > 0 && wh > 0)
       c = new TCanvas(TString::Format("c_bch1d_%d",cindex), TString::Format("c_bch1d_%d",cindex), ww, wh);
    else
       c = new TCanvas(TString::Format("c_bch1d_%d",cindex));
 
    // add c to list of objects
    fROOTObjects.push_back(c);
 
    c->cd();
 
    // set log axis
    if (flag_logx) {
       c->SetLogx();
    }
 
    // set log axis
    if (flag_logy) {
       c->SetLogy();
    }
 
    Draw(options, intervals);
 
    if (flag_rescale) {
       double top = gPad->GetTopMargin();
       double bottom = gPad->GetBottomMargin();
       double left = gPad->GetLeftMargin();
       double right = gPad->GetRightMargin();
 
       double dx = 1.-right - left;
       double dy = 1.-top-bottom;
       double ratio = dy/dx;
       double ynew = c->GetWindowWidth()/ratio;
       c->SetCanvasSize(c->GetWindowWidth(), (int) (ynew + 0.5));
       gPad->RedrawAxis();
 
       c->Modified();
       c->Update();
    }
 
    // print to file.
    c->Print(file);
 }

void BCH1D::Print	(	const char *	filename,
		std::string	options,
		double	interval,
		int	ww = `0`,
		int	wh = `0`
	)

Print distribution into a PostScript file.

Parameters

filename	Output filename
option	the draw options,

See Also: Print(const char * filename, std::string options="BTsiB3CS1D0Lmeanmode", std::vector<double> intervals=std::vector<double>(0), int ww=0, int wh=0)

Parameters

interval	an upper or lower limit
ww	canvas size in pixels along X
ww	canvas size in pixels along Y

See Also: Print(const char * filename, std::string options="BTsiB3CS1D0Lmeanmode", std::vector<double> intervals=std::vector<double>(0), int ww=0, int wh=0)

Definition at line 223 of file BCH1D.cxx.

 {
    std::vector<double> tempvec;
    tempvec.push_back(interval);
    Print(filename, options, tempvec, ww, wh);
 }

void BCH1D::SetColorScheme ( int scheme)

Sets the color scheme.

Parameters

scheme the scheme index
0 : black and white 1 : yellow-green-red 2 : blueish colors 2 : redish colors 2 : blueish colors

Definition at line 107 of file BCH1D.cxx.

 {
    fColors.clear();
 
    // color numbering
    // 0,1,2 : intervals
    // 3 : quantile lines
    // 4 : mean, mode, median
 
    if (scheme == 0) {
       fColors.push_back(18);
       fColors.push_back(16);
       fColors.push_back(14);
       fColors.push_back(kBlack);
       fColors.push_back(kBlack);
    }
    else if (scheme == 1) {
       fColors.push_back(kGreen);
       fColors.push_back(kYellow);
       fColors.push_back(kRed);
       fColors.push_back(kBlack);
       fColors.push_back(kBlack);
    }
    else if (scheme == 2) {
       fColors.push_back(kBlue+4);
       fColors.push_back(kBlue+2);
       fColors.push_back(kBlue);
       fColors.push_back(kOrange);
       fColors.push_back(kOrange);
    }
    else if (scheme == 3) {
       fColors.push_back(kRed+4);
       fColors.push_back(kRed+2);
       fColors.push_back(kRed);
       fColors.push_back(kGreen);
       fColors.push_back(kGreen);
    }
    else {
       SetColorScheme(1);
    }
 }

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::OutWarning(
          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 159 of file BCH1D.h.

160 { fMode=mode;

161 fModeFlag=1; };

void BCH1D::SetHistogram ( TH1D * hist)

inline

Sets the histogram.

Definition at line 149 of file BCH1D.h.

150 { fHistogram = hist; };

Member Data Documentation

std::vector<int> BCH1D::fColors

private

The colors of the color scheme.

Definition at line 307 of file BCH1D.h.

double BCH1D::fDefaultCLLimit

private

Default confidence level limit

Definition at line 295 of file BCH1D.h.

unsigned int BCH1D::fHCounter =0

staticprivate

Helper variable to get an unique number to be used in histogram naming

Definition at line 320 of file BCH1D.h.

TH1D* BCH1D::fHistogram

private

The 1D histogram

Definition at line 291 of file BCH1D.h.

double BCH1D::fMode

private

Global mode

Definition at line 299 of file BCH1D.h.

int BCH1D::fModeFlag

private

"Is there a global mode?" flag

Definition at line 303 of file BCH1D.h.

std::vector<TObject*> BCH1D::fROOTObjects

mutableprivate

Storage for plot objects.

Definition at line 311 of file BCH1D.h.

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

Public Member Functions

Static Private Member Functions

Private Attributes

Static Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation