BCHistogramFitter Class Reference

A class for fitting histograms with functions. More...

#include <BCHistogramFitter.h>

Inherits BCModel.

Collaboration diagram for BCHistogramFitter:

[legend]

List of all members.

Public Member Functions
Constructors and destructors
	BCHistogramFitter (TH1D *hist, TF1 *func)
	BCHistogramFitter ()
	~BCHistogramFitter ()
Member functions (miscellaneous methods)
int	CalculatePValueFast (std::vector< double > par, double &pvalue)
void	DrawFit (const char *options="", bool flaglegend=false)
int	Fit (TH1D *hist, TF1 *func)
int	Fit ()
double	FitFunction (std::vector< double > x, std::vector< double > parameters)
virtual double	LogAPrioriProbability (std::vector< double > parameters)
virtual double	LogLikelihood (std::vector< double > parameters)
Member functions (get)
TGraph *	GetErrorBand ()
TF1 *	GetFitFunction ()
TGraph *	GetGraphFitFunction ()
TH1D *	GetHistogram ()
Member functions (set)
int	SetFitFunction (TF1 *func)
void	SetFlagIntegration (bool flag)
int	SetHistogram (TH1D *hist)
Private Attributes
TGraph *	fErrorBand
TF1 *	fFitFunction
bool	fFlagIntegration
TGraph *	fGraphFitFunction
TH1D *	fHistogram

---

Detailed Description

A class for fitting histograms with functions.

Author:

Daniel Kollar

Kevin Kröninger

Version:

1.0

Date:

11.2008 This class allows fitting of a TH1D histogram using a TF1 function.

Definition at line 34 of file BCHistogramFitter.h.

---

Constructor & Destructor Documentation

BCHistogramFitter::BCHistogramFitter

(

)

The default constructor.

Definition at line 27 of file BCHistogramFitter.cxx.

                                     : BCModel("HistogramFitter")
{
   fHistogram = 0;
   fFitFunction = 0;

   // set default options and values
   this -> MCMCSetNIterationsRun(2000);
   this -> SetFillErrorBand();
   fFlagIntegration = true;
}

BCHistogramFitter::BCHistogramFitter	(	TH1D *	hist,
		TF1 *	func
	)

A constructor.

Parameters:

	hist	The histogram (TH1D).
	func	The fit function.

Definition at line 40 of file BCHistogramFitter.cxx.

                                                            : BCModel("HistogramFitter")
{
   fHistogram = 0;
   fFitFunction = 0;

   this -> SetHistogram(hist);
   this -> SetFitFunction(func);

   this -> MCMCSetNIterationsRun(2000);
   this -> SetFillErrorBand();
   fFlagIntegration = true;
}

BCHistogramFitter::~BCHistogramFitter

(

)

The default destructor.

Definition at line 147 of file BCHistogramFitter.cxx.

{}

---

Member Function Documentation

int BCHistogramFitter::CalculatePValueFast	(	std::vector< double >	par,
		double &	pvalue
	)

Calculate the p-value using fast-MCMC.

Parameters:

	par	A set of parameter values
	pvalue	The pvalue

Returns:

An error code

Definition at line 334 of file BCHistogramFitter.cxx.

{
   // check size of parameter vector
   if (par.size() != this -> GetNParameters())
   {
      BCLog::Out(BCLog::warning, BCLog::warning,"BCHistogramFitter::CalculatePValueFast() : Number of parameters is inconsistent.");
      return 0;
   }

   // check if histogram exists
   if (!fHistogram)
   {
      BCLog::Out(BCLog::warning, BCLog::warning,"BCHistogramFitter::CalculatePValueFast() : Histogram not defined.");
      return 0;
   }

   // define temporary variables
   int nbins = fHistogram -> GetNbinsX();

   std::vector <int> histogram;
   std::vector <double> expectation;
   histogram.assign(nbins, 0);
   expectation.assign(nbins, 0);

   double logp = 0;
   double logp_start = 0;
   int counter_pvalue = 0;

   // define starting distribution
   for (int ibin = 0; ibin < nbins; ++ibin)
   {
      // get bin boundaries
      double xmin = fHistogram -> GetBinLowEdge(ibin+1);
      double xmax = fHistogram -> GetBinLowEdge(ibin+2);

      // get the number of expected events
      double yexp = fFitFunction -> Integral(xmin, xmax);

      histogram[ibin]   = int(yexp);
      expectation[ibin] = yexp;

      // calculate p;
      logp += BCMath::LogPoisson(double(int(yexp)), yexp);
      logp_start += BCMath::LogPoisson(fHistogram -> GetBinContent(ibin+1), yexp);
   }

   int niter = 100000;

   // loop over iterations
   for (int iiter = 0; iiter < niter; ++iiter)
   {
      // loop over bins
      for (int ibin = 0; ibin < nbins; ++ibin)
      {
         // random step up or down in statistics for this bin
         double ptest = fRandom -> Rndm() - 0.5;

         // increase statistics by 1
         if (ptest > 0)
         {
            // calculate factor of probability
            double r = expectation.at(ibin) / double(histogram.at(ibin) + 1);

            // walk, or don't (this is the Metropolis part)
            if (fRandom -> Rndm() < r)
            {
               histogram[ibin] = histogram.at(ibin) + 1;
               logp += log(r);
            }
         }

         // decrease statistics by 1
         else if (ptest <= 0 && histogram[ibin] > 0)
         {
            // calculate factor of probability
            double r = double(histogram.at(ibin)) / expectation.at(ibin);

            // walk, or don't (this is the Metropolis part)
            if (fRandom -> Rndm() < r)
            {
               histogram[ibin] = histogram.at(ibin) - 1;
               logp += log(r);
            }
         }
         // debugKK
//       std::cout << " negative " << std::endl;

      } // end of looping over bins

      // increase counter
      if (logp < logp_start)
         counter_pvalue++;

   } // end of looping over iterations

   // calculate p-value
   pvalue = double(counter_pvalue) / double(niter);

   // no error
   return 1;
}

void BCHistogramFitter::DrawFit	(	const char *	options = "",
		bool	flaglegend = false
	)

Draw the fit in the current pad.

Definition at line 283 of file BCHistogramFitter.cxx.

{
   if (!fHistogram)
   {
      BCLog::Out(BCLog::warning, BCLog::warning,"BCHistogramFitter::DrawFit() : Histogram not defined.");
      return;
   }

   if (!fFitFunction)
   {
      BCLog::Out(BCLog::warning, BCLog::warning,"BCHistogramFitter::DrawFit() : Fit function not defined.");
      return;
   }

   // check wheather options contain "same"
   TString opt = options;
   opt.ToLower();

   // if not same, draw the histogram first to get the axes
   if(!opt.Contains("same"))
      fHistogram -> Draw(opt.Data());

   // draw the error band as central 68% probability interval
   fErrorBand = this -> GetErrorBandGraph(0.16, 0.84);
   fErrorBand -> Draw("f same");

   // now draw the histogram again since it was covered by the band
   fHistogram -> Draw(TString::Format("%ssame",opt.Data()));

   // draw the fit function on top
   fGraphFitFunction = this -> GetFitFunctionGraph( this->GetBestFitParameters() );
   fGraphFitFunction -> SetLineColor(kRed);
   fGraphFitFunction -> SetLineWidth(2);
   fGraphFitFunction -> Draw("l same");

   // draw legend
   if (flaglegend)
   {
      TLegend * legend = new TLegend(0.25, 0.75, 0.55, 0.95);
      legend -> SetBorderSize(0);
      legend -> SetFillColor(kWhite);
      legend -> AddEntry(fHistogram, "Data", "L");
      legend -> AddEntry(fGraphFitFunction, "Best fit", "L");
      legend -> AddEntry(fErrorBand, "Error band", "F");
      legend -> Draw();
   }

   gPad -> RedrawAxis();
}

int BCHistogramFitter::Fit	(	TH1D *	hist,
		TF1 *	func
	)

Performs the fit.

Parameters:

	hist	The histogram (TH1D).
	func	The fit function.

Returns:

An error code.

Definition at line 240 of file BCHistogramFitter.cxx.

{
   // set histogram
   if (hist)
      this -> SetHistogram(hist);
   else
   {
      BCLog::Out(BCLog::warning, BCLog::warning,"BCHistogramFitter::Fit() : Histogram not defined.");
      return 0;
   }

   // set function
   if (func)
      this -> SetFitFunction(func);
   else
   {
      BCLog::Out(BCLog::warning, BCLog::warning,"BCHistogramFitter::Fit() : Fit function not defined.");
      return 0;
   }

   // perform marginalization
   this -> MarginalizeAll();

   // maximize posterior probability, using the best-fit values close
   // to the global maximum from the MCMC
   this -> FindModeMinuit( this -> GetBestFitParameters() , -1);

   // calculate the p-value using the fast MCMC algorithm
   double pvalue;
   if (this -> CalculatePValueFast(this -> GetBestFitParameters(), pvalue))
      fPValue = pvalue;
   else
      BCLog::Out(BCLog::warning, BCLog::warning,"BCHistogramFitter::Fit() : Could not use the fast p-value evaluation.");

   // print summary to screen
   this -> PrintShortFitSummary();

   // no error
   return 1;
}

int BCHistogramFitter::Fit

(

)

[inline]

Performs the fit.

Returns:

An error code.

Definition at line 138 of file BCHistogramFitter.h.

         { return this -> Fit(fHistogram, fFitFunction); };

double BCHistogramFitter::FitFunction	(	std::vector< double >	x,
		std::vector< double >	parameters
	)			[virtual]

Plots the histogram

Parameters:

	options	Options for plotting.
	filename	Name of the file which the histogram is printed into. The following options are available: F : plots the fit function on top of the data E0 : plots the fit function and the 68% prob. uncertainty band of the fit function on top of the data E1 : plots the expectation from the fit function and the uncertainty bin-by-bin as error bars. Returns the y-value of the 1-dimensional fit function at an x and for a set of parameters.
	x	A vector with the x-value.
	parameters	A set of parameters.

Reimplemented from BCIntegrate.

Definition at line 230 of file BCHistogramFitter.cxx.

{
   // set the parameters of the function
   fFitFunction -> SetParameters(&params[0]);

   return fFitFunction -> Eval(x[0]) * fHistogram -> GetBinWidth( fHistogram -> FindBin(x[0]) );
}

TGraph* BCHistogramFitter::GetErrorBand

(

)

[inline]

Returns:

pointer to the error band

Definition at line 72 of file BCHistogramFitter.h.

         { return fErrorBand; };

TF1* BCHistogramFitter::GetFitFunction

(

)

[inline]

Returns:

The fit function

Definition at line 67 of file BCHistogramFitter.h.

         { return fFitFunction; };

TGraph* BCHistogramFitter::GetGraphFitFunction

(

)

[inline]

Returns:

pointer to a graph for the fit function

Definition at line 77 of file BCHistogramFitter.h.

         { return fGraphFitFunction; };

TH1D* BCHistogramFitter::GetHistogram

(

)

[inline]

Returns:

The histogram

Definition at line 62 of file BCHistogramFitter.h.

         { return fHistogram; };

double BCHistogramFitter::LogAPrioriProbability

(

std::vector< double >

parameters

)

[virtual]

The log of the prior probability. Overloaded from BCModel.

Parameters:

parameters

A vector of doubles containing the parameter values.

Reimplemented from BCModel.

Definition at line 152 of file BCHistogramFitter.cxx.

{
   // using flat probability in all parameters
   double logprob = 0.;
   for(unsigned int i=0; i < this -> GetNParameters(); i++)
      logprob -= log(this -> GetParameter(i) -> GetRangeWidth());

   return logprob;
}

double BCHistogramFitter::LogLikelihood

(

std::vector< double >

parameters

)

[virtual]

The log of the conditional probability. Overloaded from BCModel.

Parameters:

parameters

A vector of doubles containing the parameter values.

Reimplemented from BCModel.

Definition at line 164 of file BCHistogramFitter.cxx.

{
   // initialize probability
   double loglikelihood = 0;

   // set the parameters of the function
   fFitFunction -> SetParameters(&params[0]);

   // get the number of bins
   int nbins = fHistogram -> GetNbinsX();

   // get bin width
   double dx = fHistogram -> GetBinWidth(1);

   // get function value of lower bin edge
   double fedgelow = fFitFunction -> Eval(fHistogram -> GetBinLowEdge(1));
   
   // loop over all bins
   for (int ibin = 1; ibin <= nbins; ++ibin)
   {
      // get upper bin edge
      double xedgehi = fHistogram -> GetBinLowEdge(ibin) + dx;

      // get function value at upper bin edge
      double fedgehi = fFitFunction -> Eval(xedgehi);

      // get the number of observed events
      double y = fHistogram -> GetBinContent(ibin);

      double yexp = 0.;

      // use ROOT's TH1D::Integral method
      if (fFlagIntegration)
         yexp = fFitFunction -> Integral(xedgehi-dx, xedgehi);

      // use linear interpolation
      else
      {
         yexp = fedgelow * dx + (fedgehi - fedgelow) * dx / 2.;

         // make the upper edge the lower edge for the next iteration
         fedgelow = fedgehi;
      }

      // get the value of the Poisson distribution
      loglikelihood += BCMath::LogPoisson(y, yexp);
   }

// // get bin boundaries
// double xmin = fHistogram -> GetBinLowEdge(ibin);
// double xmax = fHistogram -> GetBinLowEdge(ibin+1);

// // get the number of observed events
// double y = fHistogram -> GetBinContent(ibin);

// // get the number of expected events
// double yexp = fFitFunction -> Integral(xmin, xmax);

// // get the value of the Poisson distribution
// loglikelihood += BCMath::LogPoisson(y, yexp);

   return loglikelihood;
}

int BCHistogramFitter::SetFitFunction

(

TF1 *

func

)

Parameters:

func

The fit function @ return An error code (1:pass, 0:fail).

Definition at line 109 of file BCHistogramFitter.cxx.

{
   // check if function exists
   if(!func)
   {
      BCLog::Out(BCLog::error,BCLog::error,"BCHistogramFitter::SetFitFunction() : TF1 not created.");
      return 0;
   }

   // set the function
   fFitFunction = func;

   // update the model name to contain the function name
   this -> SetName(TString::Format("HistogramFitter with %s",fFitFunction->GetName()));

   // reset parameters
   fParameterSet -> clear();

   // get the new number of parameters
   int n = func -> GetNpar();

   // add parameters
   for (int i = 0; i < n; ++i)
   {
      double xmin;
      double xmax;

      func -> GetParLimits(i, xmin, xmax);

      this -> AddParameter(func->GetParName(i), xmin, xmax);
   }

   // no error
   return 1;
}

void BCHistogramFitter::SetFlagIntegration

(

bool

flag

)

[inline]

Sets the flag for integration.
true: use ROOT's TH1D::Integrate()
false: use linear interpolation

Definition at line 101 of file BCHistogramFitter.h.

         { fFlagIntegration = flag; };

int BCHistogramFitter::SetHistogram

(

TH1D *

hist

)

Parameters:

hist

The histogram @ return An error code (1:pass, 0:fail).

Definition at line 55 of file BCHistogramFitter.cxx.

{
   // check if histogram exists
   if (!hist)
   {
      BCLog::Out(BCLog::error,BCLog::error,"BCHistogramFitter::SetHistogram() : TH1D not created.");
      return 0;
   }

   // set pointer to histogram
   fHistogram = hist;

   // create a data set. this is necessary in order to calculate the
   // error band. the data set contains as many data points as there
   // are bins. for now, the data points are empty.
   BCDataSet * ds = new BCDataSet();

   // create data points and add them to the data set.
   int nbins = fHistogram -> GetNbinsX();
   for (int i = 0; i < nbins; ++i)
   {
      BCDataPoint* dp = new BCDataPoint(2);
      ds -> AddDataPoint(dp);
   }

   // set the new data set.
   this -> SetDataSet(ds);

   // calculate the lower and upper edge in x.
   double xmin = hist -> GetBinLowEdge(1);
   double xmax = hist -> GetBinLowEdge(nbins+1);

   // calculate the minimum and maximum range in y.
   double histymin = hist -> GetMinimum();
   double histymax = hist -> GetMaximum();

   // calculate the minimum and maximum of the function value based on
   // the minimum and maximum value in y.
   double ymin = TMath::Max(0., histymin - 5.*sqrt(histymin));
   double ymax = histymax + 5.*sqrt(histymax);

   // set the data boundaries for x and y values.
   this -> SetDataBoundaries(0, xmin, xmax);
   this -> SetDataBoundaries(1, ymin, ymax);

   // set the indeces for fitting.
   this -> SetFitFunctionIndices(0, 1);

   // no error
   return 1;
}

---

Member Data Documentation

TGraph* BCHistogramFitter::fErrorBand [private]

Pointer to the error band (for legend)

Definition at line 179 of file BCHistogramFitter.h.

TF1* BCHistogramFitter::fFitFunction [private]

The fit function

Definition at line 170 of file BCHistogramFitter.h.

bool BCHistogramFitter::fFlagIntegration [private]

Flag for using the ROOT TH1D::Integral method (true), or linear interpolation (false)

Definition at line 175 of file BCHistogramFitter.h.

TGraph* BCHistogramFitter::fGraphFitFunction [private]

Pointer to a graph for displaying the fit function

Definition at line 183 of file BCHistogramFitter.h.

TH1D* BCHistogramFitter::fHistogram [private]

The histogram containing the data.

Definition at line 166 of file BCHistogramFitter.h.