The base class for all user-defined models. More...

#include <BCModel.h>

Inheritance diagram for BCModel:

Collaboration diagram for BCModel:

Public Member Functions
Constructors and destructors

	BCModel ()
	BCModel (const char *name)
virtual	~BCModel ()
Member functions (get)

std::string	GetName ()
int	GetIndex ()
double	GetModelAPrioriProbability ()
double	GetModelAPosterioriProbability ()
double	GetNormalization ()
BCDataSet *	GetDataSet ()
BCDataPoint *	GetDataPointLowerBoundaries ()
BCDataPoint *	GetDataPointUpperBoundaries ()
double	GetDataPointLowerBoundary (unsigned int index)
double	GetDataPointUpperBoundary (unsigned int index)
bool	GetFlagBoundaries ()
int	GetNDataPoints ()
BCDataPoint *	GetDataPoint (unsigned int index)
unsigned int	GetNDataPointsMinimum ()
unsigned int	GetNDataPointsMaximum ()
unsigned int	GetNParameters ()
BCParameter *	GetParameter (int index)
BCParameter *	GetParameter (const char *name)
BCParameterSet *	GetParameterSet ()
double	GetBestFitParameter (unsigned int index)
double	GetBestFitParameterError (unsigned int index)
std::vector< double >	GetBestFitParameters ()
std::vector< double >	GetBestFitParameterErrors ()
double	GetBestFitParameterMarginalized (unsigned int index)
std::vector< double >	GetBestFitParametersMarginalized ()
TH2D *	GetErrorBandXY ()
TH2D *	GetErrorBandXY_yellow (double level=.68, int nsmooth=0)
std::vector< double >	GetErrorBand (double level)
TGraph *	GetErrorBandGraph (double level1, double level2)
TGraph *	GetFitFunctionGraph (std::vector< double > parameters)
TGraph *	GetFitFunctionGraph ()
TGraph *	GetFitFunctionGraph (std::vector< double > parameters, double xmin, double xmax, int n=1000)
bool	GetFixedDataAxis (unsigned int index)
Member functions (set)

void	SetName (const char *name)
void	SetIndex (int index)
void	SetParameterSet (BCParameterSet *parset)
int	SetParameterRange (int index, double parmin, double parmax)
void	SetModelAPrioriProbability (double probability)
void	SetModelAPosterioriProbability (double probability)
void	SetNormalization (double norm)
void	SetDataSet (BCDataSet *dataset)
void	SetSingleDataPoint (BCDataPoint *datapoint)
void	SetSingleDataPoint (BCDataSet *dataset, unsigned int index)
void	SetNDataPointsMinimum (unsigned int minimum)
void	SetNDataPointsMaximum (unsigned int maximum)
void	SetDataBoundaries (unsigned int index, double lowerboundary, double upperboundary, bool fixed=false)
void	SetErrorBandContinuous (bool flag)
void	SetNbins (const char *parname, int nbins)
int	SetPrior (int index, TF1 *f)
int	SetPrior (const char name, TF1 f)
int	SetPriorGauss (int index, double mean, double sigma)
int	SetPriorGauss (const char *name, double mean, double sigma)
int	SetPriorGauss (int index, double mean, double sigmadown, double sigmaup)
int	SetPriorGauss (const char *name, double mean, double sigmadown, double sigmaup)
int	SetPrior (int index, TH1 *h, bool flag=false)
int	SetPrior (const char name, TH1 h, bool flag=false)
int	SetPriorConstant (int index)
int	SetPriorConstant (const char *name)
int	SetPriorConstantAll ()
Member functions (miscellaneous methods)

int	AddParameter (const char *name, double lowerlimit, double upperlimit)
int	AddParameter (BCParameter *parameter)
double	APrioriProbability (std::vector< double > parameters)
virtual double	LogAPrioriProbability (std::vector< double > parameters)
double	Likelihood (std::vector< double > parameter)
virtual double	LogLikelihood (std::vector< double > parameter)
double	ProbabilityNN (std::vector< double > parameter)
double	LogProbabilityNN (std::vector< double > parameter)
double	Probability (std::vector< double > parameter)
double	LogProbability (std::vector< double > parameter)
double	ConditionalProbabilityEntry (BCDataPoint *datapoint, std::vector< double > parameters)
virtual double	LogConditionalProbabilityEntry (BCDataPoint *, std::vector< double >)
virtual double	SamplingFunction (std::vector< double > parameters)
double	Eval (std::vector< double > parameters)
double	LogEval (std::vector< double > parameters)
double	EvalSampling (std::vector< double > parameters)
double	Normalize ()
int	CheckParameters (std::vector< double > parameters)
void	FindMode (std::vector< double > start=std::vector< double >(0))
void	FindModeMinuit (std::vector< double > start=std::vector< double >(0), int printlevel=1)
void	WriteMode (const char *file)
int	ReadMode (const char *file)
int	ReadMarginalizedFromFile (const char *file)
int	ReadErrorBandFromFile (const char *file)
int	MarginalizeAll ()
BCH1D *	GetMarginalized (BCParameter *parameter)
BCH1D *	GetMarginalized (const char *name)
BCH2D *	GetMarginalized (BCParameter parameter1, BCParameter parameter2)
BCH2D *	GetMarginalized (const char name1, const char name2)
int	PrintAllMarginalized1D (const char *filebase)
int	PrintAllMarginalized2D (const char *filebase)
int	PrintAllMarginalized (const char *file, unsigned int hdiv=1, unsigned int ndiv=1)
virtual void	CorrelateDataPointValues (std::vector< double > &x)
double	GetPvalueFromChi2 (std::vector< double > par, int sigma_index)
double	GetPvalueFromChi2Johnson (std::vector< double > par)
double	GetPvalueFromKolmogorov (const std::vector< double > &par, int index)
double	GetChi2Johnson (std::vector< double > par, const int nBins=-1)
double	GetAvalueFromChi2Johnson (TTree tree, TH1D distribution=0)
double	GetPvalueFromChi2NDoF (std::vector< double > par, int sigma_index)
BCH1D *	CalculatePValue (std::vector< double > par, bool flag_histogram=false)
double	GetPValue ()
double	GetPValueNDoF ()
double	GetChi2NDoF ()
std::vector< double >	GetChi2Runs (int dataIndex, int sigmaIndex)
void	SetGoFNIterationsMax (int n)
void	SetGoFNIterationsRun (int n)
void	SetGoFNChains (int n)
double	HessianMatrixElement (BCParameter parameter1, BCParameter parameter2, std::vector< double > point)
void	PrintSummary ()
void	PrintResults (const char *file)
void	PrintShortFitSummary (int chi2flag=0)
void	PrintHessianMatrix (std::vector< double > parameters)
void	FixDataAxis (unsigned int index, bool fixed)
virtual double	CDF (const std::vector< double > &, int, bool)
int	ResetResults ()
Protected Attributes
double	fChi2NDoF
BCDataSet *	fDataSet
int	fGoFNChains
int	fGoFNIterationsMax
int	fGoFNIterationsRun
int	fIndex
bool	flag_ConditionalProbabilityEntry
bool	flag_discrete
double	fModelAPosteriori
double	fModelAPriori
std::string	fName
unsigned int	fNDataPointsMaximum
unsigned int	fNDataPointsMinimum
BCParameterSet *	fParameterSet
bool	fPriorConstantAll
double	fPriorConstantValue
std::vector< TNamed * >	fPriorContainer
std::vector< bool >	fPriorContainerConstant
std::vector< bool >	fPriorContainerInterpolate
double	fPValue
double	fPValueNDoF
Private Member Functions
int	CompareStrings (const char string1, const char string2)
int	NumberBins ()
void	RecalculatePriorConstant ()
BCDataPoint *	VectorToDataPoint (std::vector< double > data)
Private Attributes
double	fNormalization

Detailed Description

The base class for all user-defined models.

Author:: Daniel Kollar; Kevin Kröninger

Version:: 1.0

Date:: 08.2008 This class represents a model. It contains a container of parameters, their prior distributions and the conditional probabilities given those parameters. The methods which implement the prior and conditional probabilities have to be overloaded by the user in the user defined model class which will inherit from this class.

Definition at line 53 of file BCModel.h.

Constructor & Destructor Documentation

BCModel::BCModel ( )

The default constructor.

Definition at line 78 of file BCModel.cxx.

   : BCIntegrate()
{
   fNormalization = -1.;
   fDataSet = 0;
   fParameterSet = new BCParameterSet();

   fIndex = -1;
   fPValue = -1;
   fPValueNDoF = -1;
   fChi2NDoF = -1;

   fName = "model";
   fDataPointUpperBoundaries = 0;
   fDataPointLowerBoundaries = 0;

   flag_ConditionalProbabilityEntry = true;

   fGoFNChains = 5;
   fGoFNIterationsMax = 100000;
   fGoFNIterationsRun = 2000;

   flag_discrete = false;

   fPriorConstantAll = false;
   fPriorConstantValue = 0;
}

BCModel::BCModel ( const char * name )

A constructor.

Parameters:

name

The name of the model

Definition at line 47 of file BCModel.cxx.

   : BCIntegrate()
{
   fNormalization = -1.;
   fDataSet = 0;
   fParameterSet = new BCParameterSet;

   fIndex = -1;
   fPValue = -1;
   fPValueNDoF = -1;
   fChi2NDoF = -1;

   fName = (char *) name;
   flag_ConditionalProbabilityEntry = true;

   fDataPointUpperBoundaries = 0;
   fDataPointLowerBoundaries = 0;

   fErrorBandXY = 0;

   fGoFNChains = 5;
   fGoFNIterationsMax = 100000;
   fGoFNIterationsRun = 2000;

   flag_discrete = false;

   fPriorConstantAll = false;
   fPriorConstantValue = 0;
}

BCModel::~BCModel ( ) [virtual]

The default destructor.

Definition at line 107 of file BCModel.cxx.

{
   for (unsigned int i = 0; i < GetNParameters(); ++i) 
      delete fPriorContainer[i]; 
   fPriorContainer.clear();

   delete fParameterSet;

   delete fDataPointLowerBoundaries;
   delete fDataPointUpperBoundaries;
}

Member Function Documentation

int BCModel::AddParameter	(	const char *	name,
		double	lowerlimit,
		double	upperlimit
	)

Adds a parameter to the parameter set

Parameters:

	name	The name of the parameter
	lowerlimit	The lower limit of the parameter values
	upperlimit	The upper limit of the parameter values

See also:: AddParameter(BCParameter* parameter);

Definition at line 469 of file BCModel.cxx.

{
   // create new parameter
   BCParameter * parameter = new BCParameter(name, lowerlimit, upperlimit);

   int flag_ok = AddParameter(parameter);
   if (flag_ok)
      delete parameter;

   return flag_ok;
}

int BCModel::AddParameter ( BCParameter * parameter )

Adds a parameter to the model.

Parameters:

parameter

A model parameter

See also:: AddParameter(const char * name, double lowerlimit, double upperlimit);

Definition at line 482 of file BCModel.cxx.

{
   // check if parameter set exists
   if (!fParameterSet) {
      BCLog::OutError("BCModel::AddParameter : Parameter set does not exist");
      return ERROR_PARAMETERSETDOESNOTEXIST;
   }

   // check if parameter with same name exists
   int flag_exists = 0;
   for (unsigned int i = 0; i < GetNParameters(); i++)
      if (CompareStrings(parameter->GetName().data(), GetParameter(i)->GetName().data()) == 0)
         flag_exists = -1;

   if (flag_exists < 0) {
      BCLog::OutError(Form(
            "BCModel::AddParameter : Parameter with name %s exists already. ",
            parameter->GetName().data()));
      return ERROR_PARAMETEREXISTSALREADY;
   }

   // define index of new parameter
   unsigned int index = fParameterSet->size();
   parameter->SetIndex(index);

   // add parameter to parameter container
   fParameterSet->push_back(parameter);

   // add empty object to prior container
   fPriorContainer.push_back(0); 

   // don't interpolate the prior histogram by default
   fPriorContainerInterpolate.push_back(false);

   // prior assumed to be non-constant in general case
   fPriorContainerConstant.push_back(false);

   // add parameters to integration methods
   SetParameters(fParameterSet);

   // reset results
   ResetResults();

   return 0;
}

double BCModel::APrioriProbability ( std::vector< double > parameters ) [inline]

Returns the prior probability.

Parameters:

parameters

A set of parameter values

Returns:: The prior probability p(parameters)

See also:: GetPrior(std::vector <double> parameters)

Definition at line 444 of file BCModel.h.

         { return exp( this->LogAPrioriProbability(parameters) ); };

BCH1D * BCModel::CalculatePValue	(	std::vector< double >	par,
		bool	flag_histogram = `false`
	)

Definition at line 1576 of file BCModel.cxx.

{
   BCH1D * hist = 0;

   // print log
   BCLog::OutSummary("Do goodness-of-fit-test");

   // create model test
   BCGoFTest * goftest = new BCGoFTest("modeltest");

   // set this model as the model to be tested
   goftest->SetTestModel(this);

   // set the point in parameter space which is tested an initialize
   // the model testing
   if (!goftest->SetTestPoint(par))
      return 0;

   // disable the creation of histograms to save _a lot_ of memory
   // (histograms are not needed for p-value calculation)
   goftest->MCMCSetFlagFillHistograms(false);

   // set parameters of the MCMC for the GoFTest
   goftest->MCMCSetNChains(fGoFNChains);
   goftest->MCMCSetNIterationsMax(fGoFNIterationsMax);
   goftest->MCMCSetNIterationsRun(fGoFNIterationsRun);

   // get p-value
   fPValue = goftest->GetCalculatedPValue(flag_histogram);

   // get histogram
   if (flag_histogram) {
      hist = new BCH1D();
      hist->SetHistogram(goftest->GetHistogramLogProb());
   }

   // delete model test
   delete goftest;

   // return histogram
   return hist;
}

virtual double BCModel::CDF	(	const std::vector< double > &	,
		int	,
		bool
	)			`[inline, virtual]`

1dim cumulative distribution function of the probability to get the data f(x_i|param) for a single measurement, assumed to be of identical functional form for all measurements

Parameters:

	parameters	The parameter values at which point to compute the cdf
	index	The data point index starting at 0,1...N-1
	lower	only needed for discrete distributions! Return the CDF for the count one less than actually observed, e.g. in Poisson process, if 3 actually observed, then CDF(2) is returned

Reimplemented in BCGraphFitter, and BCHistogramFitter.

Definition at line 751 of file BCModel.h.

      {return 0.0;}

int BCModel::CheckParameters ( std::vector< double > parameters )

Checks if a set of parameters values is within the given range.

Parameters:

parameters

A set of parameter values

Returns:: Error code (0: OK, -1 length of parameters not correct, -2 values not within range)

Definition at line 663 of file BCModel.cxx.

{
   // check if vectors are of equal size
   if (!parameters.size() == fParameterSet->size())
      return ERROR_INVALIDNUMBEROFPARAMETERS;

   // check if parameters are within limits
   for (unsigned int i = 0; i < fParameterSet->size(); i++) {
      BCParameter * modelparameter = fParameterSet->at(i);

      if (modelparameter->GetLowerLimit() > parameters.at(i)
          || modelparameter->GetUpperLimit() < parameters.at(i)) {
         BCLog::OutError(Form(
               "BCModel::CheckParameters : Parameter %s not within limits.",
               fParameterSet-> at(i)-> GetName().data()));
         return ERROR_PARAMETERNOTWITHINRANGE;
      }
   }

   return 0;
}

int BCModel::CompareStrings	(	const char *	string1,
		const char *	string2
	)			`[private]`

Compares to strings

Definition at line 2286 of file BCModel.cxx.

{
   int flag_same = 0;

   if (strlen(string1) != strlen(string2))
      return -1;

   for (int i = 0; i < int(strlen(string1)); i++)
      if (string1[i] != string2[i])
         flag_same = -1;

   return flag_same;
}

double BCModel::ConditionalProbabilityEntry	(	BCDataPoint *	datapoint,
		std::vector< double >	parameters
	)			`[inline]`

Returns a conditional probability. Method needs to be overloaded by the user.

Parameters:

	datapoint	A data point
	parameters	A set of parameter values

Returns:: The conditional probability p(datapoint|parameters)

See also:: GetConditionalEntry(BCDataPoint* datapoint, std::vector <double> parameters)

Definition at line 503 of file BCModel.h.

         { return exp( this->LogConditionalProbabilityEntry(datapoint, parameters) ); };

void BCModel::CorrelateDataPointValues ( std::vector< double > & x ) [virtual]

Constrains a data point

Parameters:

x

A vector of double

Definition at line 1620 of file BCModel.cxx.

{
   // ...
}

double BCModel::Eval ( std::vector< double > parameters ) [inline, virtual]

Overloaded function to evaluate integral.

Reimplemented from BCIntegrate.

Definition at line 525 of file BCModel.h.

         { return exp( this->LogEval(parameters) ); };

double BCModel::EvalSampling ( std::vector< double > parameters ) [virtual]

Overloaded function to evaluate integral.

Reimplemented from BCIntegrate.

Definition at line 621 of file BCModel.cxx.

{
   return SamplingFunction(parameters);
}

void BCModel::FindMode ( std::vector< double > start = std::vector<double>(0) )

Do the mode finding using a method set via SetOptimizationMethod. Default is Minuit. The mode can be extracted using the GetBestFitParameters() method.

A starting point for the mode finding can be specified for Minuit. If not specified, Minuit default will be used (center of the parameter space).

If running mode finding after the MCMC it is a good idea to specify the mode obtained from MCMC as a starting point for the Minuit minimization. MCMC will find the location of the global mode and Minuit will converge to the mode precisely. The commands are:

			model -> MarginalizeAll();
			model -> FindMode( model -> GetBestFitParameters() );

startinf point of Minuit minimization

Definition at line 686 of file BCModel.cxx.

{
   if(fParameterSet->size()<1) {
      BCLog::OutError(Form("FindMode : No parameters defined in model \'%s\'. Aborting.",GetName().data()));
      return;
   }

   // this implementation is CLEARLY not good we have to work on this.

   BCLog::OutSummary(Form("Model \'%s\': Finding mode using %s",GetName().data(), DumpOptimizationMethod().c_str()));

   // synchronize parameters in BCIntegrate
   SetParameters(fParameterSet);

   switch (GetOptimizationMethod()) {
      case BCIntegrate::kOptSA:
         FindModeSA(start);
         return;

      case BCIntegrate::kOptMinuit: {
         int printlevel = -1;
         if (BCLog::GetLogLevelScreen() <= BCLog::detail)
            printlevel = 0;
         if (BCLog::GetLogLevelScreen() <= BCLog::debug)
            printlevel = 1;
         BCIntegrate::FindModeMinuit(start, printlevel);
         return;
      }

      case BCIntegrate::kOptMetropolis:
          MarginalizeAll();
         return;

      default:
         BCLog::OutError(Form("BCModel::FindMode : Invalid mode finding method: %d",GetOptimizationMethod()));
   }


   return;
}

void BCModel::FindModeMinuit	(	std::vector< double >	start = `std::vector<double>(0)`,
		int	printlevel = `1`
	)			`[virtual]`

Does the mode finding using Minuit. If starting point is not specified, finding will start from the center of the parameter space.

Parameters:

	start	point in parameter space from which the mode finding is started.
	printlevel	The print level.

Reimplemented from BCIntegrate.

Definition at line 728 of file BCModel.cxx.

{
   if(fParameterSet->size()<1) {
      BCLog::OutError(Form("FindModeMinuit : No parameters defined in model \'%s\'. Aborting.",GetName().data()));
      return;
   }

   // synchronize parameters in BCIntegrate
   SetParameters(fParameterSet);
   BCIntegrate::FindModeMinuit(start, printlevel);
}

void BCModel::FixDataAxis	(	unsigned int	index,
		bool	fixed
	)

Definition at line 1671 of file BCModel.cxx.

{
   // check if index is within range
   if (index < 0 || index > fDataSet->GetDataPoint(0)->GetNValues()) {
      BCLog::OutError("BCModel::FixDataAxis : Index out of range.");
      return;
   }

   if (fDataFixedValues.size() == 0)
      fDataFixedValues.assign(fDataSet->GetDataPoint(0)->GetNValues(),
            false);

   fDataFixedValues[index] = fixed;
}

double BCModel::GetAvalueFromChi2Johnson	(	TTree *	tree,
		TH1D *	distribution = `0`
	)

Calculate the A-value, a summary statistic. It computes the frequency that a Chi2 value determined from the data by Johnson's binning prescription is larger than a value sampled from the reference chi2 distribution. They out from one chain is used. A=1/2 provides no evidence against the null hypothesis. Compare Johnson, V.E. A Bayesian chi2 Test for Goodness-of-Fit. The Annals of Statistics 32, 2361-2384(2004).

Parameters:

	par	tree contains the samples of posterior of the parameters
	par	histogram filled by function with distribution of p values

Definition at line 1440 of file BCModel.cxx.

{
   // model parameters
   int nPar = (int)GetNParameters();
   std::vector<double> param(nPar);

   //parameters saved in branches should be the same
   int nParBranches = -1;
   tree->SetBranchAddress("fNParameters", &nParBranches);

   //assume all events have same number of parameters, so check only first
   tree->GetEntry(0);
   if (nParBranches != nPar) {
      BCLog::OutError(Form(
            "Cannot compute A value: number of parameters in tree (%d)"
               "doesn't match  number of parameters in model (%d)",
            nParBranches, nPar));
      return -1.;
   }

   // buffer to construct correct branchnames for parameters, e.g. "fParameter3"
   char * branchName = new char[10 + nPar];
   // set up variables filled for each sample of parameters
   // assume same order as in model
   for (int i = 0; i < (int) nPar; i++) {
      +sprintf(branchName, "fParameter%d", i);
      tree->SetBranchAddress(branchName, &param[i]);
   }

   // get the p value from Johson's definition for each param from posterior
   long nEntries = tree->GetEntries();

   // RN ~ chi2 with K-1 DoF needed for comparison
   std::vector<double> randoms(nEntries);
   int K = NumberBins();
   BCMath::RandomChi2(randoms, K - 1);

   // number of Johnson chi2 values bigger than reference
   int nBigger = 0;
   for (int i = 0; i < nEntries; i++) {
      tree->GetEntry(i);
      double chi2 = GetChi2Johnson(param);
      if (distribution != 0)
         distribution->Fill(chi2);
      // compare to set of chi2 variables
      if (chi2 > randoms[i])
         nBigger++;
   }

   return nBigger / double(nEntries);
}

double BCModel::GetBestFitParameter ( unsigned int index )

Returns the value of a parameter (defined by index) at the global mode of the posterior pdf.

Parameters:

index

index of the parameter.

Returns:: best fit value of the parameter or -1e+111 on error or center of the range if mode finding not yer run

Definition at line 180 of file BCModel.cxx.

{
   if(index >= GetNParameters()) {
      BCLog::OutError("BCModel::GetBestFitParameter : Parameter index out of range, returning -1e+111.");
      return -1e+111;
   }

   if(fBestFitParameters.size()==0) {
      BCLog::OutError("BCModel::GetBestFitParameter : Mode finding not yet run, returning center of the range.");
      return (GetParameter(index)->GetUpperLimit() + GetParameter(index)->GetLowerLimit() ) / 2.;
   }

   return fBestFitParameters[index];
}

double BCModel::GetBestFitParameterError ( unsigned int index )

Returns the error on the value of a parameter (defined by index) at the global mode of the posterior pdf.

Parameters:

index

index of the parameter.

Returns:: error on the best fit value of the parameter or -1 if undefined

Definition at line 196 of file BCModel.cxx.

{
   if(index >= GetNParameters()) {
      BCLog::OutError("BCModel::GetBestFitParameterError : Parameter index out of range, returning -1.");
      return -1;
   }

   if(fBestFitParameterErrors.size()==0) {
      BCLog::OutError("BCModel::GetBestFitParameterError : Mode finding not yet run, returning -1.");
      return -1.;
   }

   if(fBestFitParameterErrors[index]<0.)
      BCLog::OutWarning("BCModel::GetBestFitParameterError : Parameter error not available, returning -1.");

   return fBestFitParameterErrors[index];
}

std::vector<double> BCModel::GetBestFitParameterErrors ( ) [inline]

Definition at line 196 of file BCModel.h.

         { return fBestFitParameterErrors; };

double BCModel::GetBestFitParameterMarginalized ( unsigned int index )

Returns the value of a particular parameter (defined by index) at the modes of the marginalized posterior pdfs.

Parameters:

index

index of the parameter.

Returns:: best fit parameter or -1e+111 on error or center of the range if marginalization not yet run

Definition at line 215 of file BCModel.cxx.

{
   if(index >= GetNParameters()) {
      BCLog::OutError("BCModel::GetBestFitParameterMarginalized : Parameter index out of range, returning -1e+111.");
      return -1e+111;
   }

   if(fBestFitParametersMarginalized.size()==0) {
      BCLog::OutError("BCModel::GetBestFitParameterMarginalized : MCMC not yet run, returning center of the range.");
      return (GetParameter(index)->GetUpperLimit() + GetParameter(index)->GetLowerLimit() ) / 2.;
   }

   return fBestFitParametersMarginalized[index];
}

std::vector<double> BCModel::GetBestFitParameters ( ) [inline]

Returns the set of values of the parameters at the global mode of the posterior pdf.

Returns:: The best fit parameters

Definition at line 193 of file BCModel.h.

         { return fBestFitParameters; };

std::vector<double> BCModel::GetBestFitParametersMarginalized ( ) [inline]

Returns the set of values of the parameters at the modes of the marginalized posterior pdfs.

Returns:: best fit parameters

Definition at line 210 of file BCModel.h.

         { return fBestFitParametersMarginalized; };

double BCModel::GetChi2Johnson	(	std::vector< double >	par,
		const int	nBins = `-1`
	)

Calculate Chi2 (also called R^{B}) for arbitrary problems with binned data using the definition (3) from Johnson, V.E. A Bayesian chi2 Test for Goodness-of-Fit. The Annals of Statistics 32, 2361-2384(2004).

Parameters:

	par	Parameter set for the calculation of the likelihood
	nBins	how many bins to use for the data, for negative an adapted rule \ of thumb by Wald(1942) is used, with at least three bins

Definition at line 1327 of file BCModel.cxx.

{
   typedef unsigned int uint;

   // number of observations
   int n = GetNDataPoints();

   if (nBins < 0)
      nBins = NumberBins();

   // fixed width quantiles, including final point!
   std::vector<double> a;
   for (int i = 0; i <= nBins; i++)
      a.push_back(i / double(nBins));

   // determine the bin counts and fill the histogram with data using the CDF
   TH1I * hist = new TH1I("h1", "h1 title", nBins, 0., 1.);

   // discrete case requires randomization to allocate counts of bins that cover more
   // than one quantile
   if (flag_discrete) {
      // loop over observations, each may have different likelihood and CDF
      for (int j = 0; j < n; j++) {
         // actual value
            double CDFval = CDF(par, j, false);
         // for the bin just before
         double CDFlower = CDF(par, j, true);

         // what quantiles q are covered, count from q_1 to q_{nBins}
         int qMax = 1;
         for (int i = 0; i < nBins; i++) {
            if (CDFval > a[i])
               qMax = i + 1;
            else
               break;
         }
         int qMin = 1;
         for (int i = 0; i < nBins; i++) {
            if (CDFlower > a[i])
               qMin = i + 1;
            else
               break;
         }

         // simplest case: observation bin entirely contained in one quantile
         if (qMin == qMax) {
            // watch out for overflow because CDF exactly 1
            if (CDFval < 1)
               hist->Fill(CDFval);
            else
               hist->AddBinContent(qMax);

            continue; // this observation finished
         }

         // if more than quantile is covered need more work:
         // determine probabilities of this observation to go for each quantile covered
         // as follows: If each quantile has size 0.25 and the CDF(integral of likelihood)
         // for current observation gives gives 0.27, but for observation-1 we would have
         // 0.20, then 5/7 of the 7% go for first quantile and 2/7 for the second.
         // This extend to bins covering more than two quantiles
         std::vector<double> prob;
         // normalization
         double norm = 1 / double(CDFval - CDFlower);

         for (int i = 0; i < (qMax - qMin + 1); i++) {
            if (i == 0) {
               prob.push_back(norm * (a[qMin] - CDFlower));
               continue;
            }
            if (i == (qMax - qMin)) {
               prob.push_back(norm * (CDFval - a[qMax - 1]));
               continue;
            }
            // default case
            prob.push_back(norm * (a[i] - a[i - 1]));
         }
         // have distribution, use inverse-transform method
         double U = fRandom->Rndm();
         // build up the integral (CDF)
         for (uint i = 1; i < prob.size(); i++)
            prob[i] += prob[i - 1];
         // and search with linear comput. complexity
         for (uint i = 0; i < prob.size(); i++) {
            // we finally allocate the count, as center of quantile
            if (U <= prob[i]) {
               hist->Fill((a[qMin + i - 1] + a[qMin + i]) / 2.);
               break;
            }
         }
      }
   }
   else { //continuous case is simple
      for (int j = 0; j < n; j++)
            hist->Fill(CDF(par, j, false));
   }

   // calculate chi^2
   double chi2 = 0.;
   double mk, pk;
   double N = double(n);
   for (int i = 1; i <= nBins; i++) {
      mk = hist->GetBinContent(i);
      pk = a[i] - a[i - 1];
      chi2 += (mk - N * pk) * (mk - N * pk) / (N * pk);
   }

   delete hist;

   return chi2;
}

double BCModel::GetChi2NDoF ( ) [inline]

Definition at line 687 of file BCModel.h.

         { return fChi2NDoF; };

std::vector< double > BCModel::GetChi2Runs	(	int	dataIndex,
		int	sigmaIndex
	)

For a Gaussian problem, calculate the chi2 of the longest run of consecutive values above/below the expected values

Parameters:

	dataIndex	component of datapoint with the observed value
	sigmaIndex	component of datapoint with uncertainty

Definition at line 1311 of file BCModel.cxx.

{
   std::vector<double> x;
   return x;
}

BCDataPoint * BCModel::GetDataPoint ( unsigned int index )

Parameters:

index

The index of the data point.

Returns:: The data point in the current data set at index

Definition at line 134 of file BCModel.cxx.

{
   if (fDataSet)
      return fDataSet->GetDataPoint(index);

   BCLog::OutWarning("BCModel::GetDataPoint : No data set defined.");
   return 0;
}

BCDataPoint* BCModel::GetDataPointLowerBoundaries ( ) [inline]

Returns:: The lower boundaries of possible data values.

Definition at line 111 of file BCModel.h.

         { return fDataPointLowerBoundaries; };

double BCModel::GetDataPointLowerBoundary ( unsigned int index ) [inline]

Parameters:

index

The index of the variable.

Returns:: The lower boundary of possible data values for a particular variable.

Definition at line 122 of file BCModel.h.

         { return fDataPointLowerBoundaries -> GetValue(index); };

BCDataPoint* BCModel::GetDataPointUpperBoundaries ( ) [inline]

Returns:: The upper boundaries of possible data values.

Definition at line 116 of file BCModel.h.

         { return fDataPointUpperBoundaries; };

double BCModel::GetDataPointUpperBoundary ( unsigned int index ) [inline]

Parameters:

index

The index of the variable.

Returns:: The upper boundary of possible data values for a particular variable.

Definition at line 128 of file BCModel.h.

         { return fDataPointUpperBoundaries -> GetValue(index); };

BCDataSet* BCModel::GetDataSet ( ) [inline]

Returns:: The data set.

Definition at line 106 of file BCModel.h.

         { return fDataSet; };

std::vector< double > BCModel::GetErrorBand ( double level )

Returns a vector of y-values at a certain probability level.

Parameters:

level

The level of probability

Returns:: vector of y-values

Definition at line 244 of file BCModel.cxx.

{
   std::vector<double> errorband;

   if (!fErrorBandXY)
      return errorband;

   int nx = fErrorBandXY->GetNbinsX();
   errorband.assign(nx, 0.);

   // loop over x and y bins
   for (int ix = 1; ix <= nx; ix++) {
      TH1D * temphist = fErrorBandXY->ProjectionY("temphist", ix, ix);

      int nprobSum = 1;
      double q[1];
      double probSum[1];
      probSum[0] = level;

      temphist->GetQuantiles(nprobSum, q, probSum);

      errorband[ix - 1] = q[0];
   }

   return errorband;
}

TGraph * BCModel::GetErrorBandGraph	(	double	level1,
		double	level2
	)

Definition at line 272 of file BCModel.cxx.

{
   if (!fErrorBandXY)
      return 0;

   // define new graph
   int nx = fErrorBandXY->GetNbinsX();

   TGraph * graph = new TGraph(2 * nx);
   graph->SetFillStyle(1001);
   graph->SetFillColor(kYellow);

   // get error bands
   std::vector<double> ymin = GetErrorBand(level1);
   std::vector<double> ymax = GetErrorBand(level2);

   for (int i = 0; i < nx; i++) {
      graph->SetPoint(i, fErrorBandXY->GetXaxis()->GetBinCenter(i + 1), ymin[i]);
      graph->SetPoint(nx + i, fErrorBandXY->GetXaxis()->GetBinCenter(nx - i), ymax[nx - i - 1]);
   }

   return graph;
}

TH2D* BCModel::GetErrorBandXY ( ) [inline]

Returns:: The 2-d histogram of the error band.

Definition at line 215 of file BCModel.h.

         { return fErrorBandXY; };

TH2D * BCModel::GetErrorBandXY_yellow	(	double	level = `.68`,
		int	nsmooth = `0`
	)

Definition at line 297 of file BCModel.cxx.

{
   if (!fErrorBandXY)
      return 0;

   int nx = fErrorBandXY->GetNbinsX();
   int ny = fErrorBandXY->GetNbinsY();

   // copy existing histogram
   TH2D * hist_tempxy = (TH2D*) fErrorBandXY->Clone(
         TString::Format("%s_sub_%f.2", fErrorBandXY->GetName(), level));
   hist_tempxy->Reset();
   hist_tempxy->SetFillColor(kYellow);

   // loop over x bins
   for (int ix = 1; ix < nx; ix++) {
      BCH1D * hist_temp = new BCH1D();

      TH1D * hproj = fErrorBandXY->ProjectionY("temphist", ix, ix);
      if (nsmooth > 0)
         hproj->Smooth(nsmooth);

      hist_temp->SetHistogram(hproj);

      TH1D * hist_temp_yellow = hist_temp->GetSmallestIntervalHistogram(level);

      for (int iy = 1; iy <= ny; ++iy)
         hist_tempxy->SetBinContent(ix, iy, hist_temp_yellow->GetBinContent(iy));

      delete hist_temp_yellow;
      delete hist_temp;
   }

   return hist_tempxy;
}

TGraph * BCModel::GetFitFunctionGraph	(	std::vector< double >	parameters,
		double	xmin,
		double	xmax,
		int	n = `1000`
	)

Definition at line 359 of file BCModel.cxx.

{
   // define new graph
   TGraph * graph = new TGraph(n + 1);

   double dx = (xmax - xmin) / (double) n;

   // loop over x values
   for (int i = 0; i <= n; i++) {
      double x = (double) i * dx + xmin;
      std::vector<double> xvec;
      xvec.push_back(x);
      double y = FitFunction(xvec, parameters);

      xvec.clear();

      graph->SetPoint(i, x, y);
   }

   return graph;
}

TGraph * BCModel::GetFitFunctionGraph ( std::vector< double > parameters )

Definition at line 334 of file BCModel.cxx.

{
   if (!fErrorBandXY)
      return 0;

   // define new graph
   int nx = fErrorBandXY->GetNbinsX();
   TGraph * graph = new TGraph(nx);

   // loop over x values
   for (int i = 0; i < nx; i++) {
      double x = fErrorBandXY->GetXaxis()->GetBinCenter(i + 1);

      std::vector<double> xvec;
      xvec.push_back(x);
      double y = FitFunction(xvec, parameters);
      xvec.clear();

      graph->SetPoint(i, x, y);
   }

   return graph;
}

TGraph* BCModel::GetFitFunctionGraph ( ) [inline]

Definition at line 230 of file BCModel.h.

         { return this -> GetFitFunctionGraph(this -> GetBestFitParameters()); };

bool BCModel::GetFixedDataAxis ( unsigned int index )

Definition at line 1687 of file BCModel.cxx.

{
   // check if index is within range
   if (index < 0 || index > fDataSet->GetDataPoint(0)->GetNValues()) {
      BCLog::OutError("BCModel::GetFixedDataAxis : Index out of range.");
      return false;
   }

   return fDataFixedValues[index];
}

bool BCModel::GetFlagBoundaries ( )

Definition at line 382 of file BCModel.cxx.

{
   if (!fDataPointLowerBoundaries)
      return false;

   if (!fDataPointUpperBoundaries)
      return false;

   if (fDataPointLowerBoundaries->GetNValues() != fDataSet->GetDataPoint(0)->GetNValues())
      return false;

   if (fDataPointUpperBoundaries->GetNValues() != fDataSet->GetDataPoint(0)->GetNValues())
      return false;

   return true;
}

int BCModel::GetIndex ( ) [inline]

Returns:: The index of the model.

Definition at line 86 of file BCModel.h.

         { return fIndex; };

BCH1D * BCModel::GetMarginalized ( BCParameter * parameter )

If MarginalizeAll method was used, the individual marginalized distributions with respect to one parameter can be retrieved using this method.

Parameters:

parameter

Model parameter

Returns:: 1D marginalized probability

Definition at line 868 of file BCModel.cxx.

{
   if (!parameter) {
// don't print any error message, should be done upstream
//      BCLog::OutError("BCModel::GetMarginalized : Parameter does not exist.");
      return 0;
   }

   int index = parameter->GetIndex();
   if (fMCMCFlagsFillHistograms[index] == false) {
// don't print any error message, should be done upstream
//      BCLog::OutError(Form(
//            "BCModel::GetMarginalized : Distribuion for '%s' not filled.",
//            parameter->GetName().data()));
      return 0;
   }

   // get histogram
   TH1D * hist = MCMCGetH1Marginalized(index);
   if (!hist)
      return 0;

   // set axis labels
   hist->SetName(Form("hist_%s_%s", GetName().data(), parameter->GetName().data()));
   hist->SetXTitle(parameter->GetName().data());
   hist->SetYTitle(Form("p(%s|data)", parameter->GetName().data()));
   hist->SetStats(kFALSE);

   // set histogram
   BCH1D * hprob = new BCH1D();
   hprob->SetHistogram(hist);

   // set best fit parameter
   double bestfit = hprob->GetMode();

   if (fBestFitParametersMarginalized.size() == 0)
      for (unsigned int i = 0; i < GetNParameters(); i++)
         fBestFitParametersMarginalized.push_back(0.);

   fBestFitParametersMarginalized[index] = bestfit;

   hprob->SetGlobalMode(fBestFitParameters.at(index));

   return hprob;
}

BCH1D* BCModel::GetMarginalized ( const char * name ) [inline]

Definition at line 603 of file BCModel.h.

         { return this -> GetMarginalized(this -> GetParameter(name)); };

BCH2D * BCModel::GetMarginalized	(	BCParameter *	parameter1,
		BCParameter *	parameter2
	)

If MarginalizeAll method was used, the individual marginalized distributions with respect to two parameters can be retrieved using this method.

Parameters:

	parameter1	First parameter
	parameter2	Second parameter

Returns:: 2D marginalized probability

Definition at line 1239 of file BCModel.cxx.

{
   int index1 = par1->GetIndex();
   int index2 = par2->GetIndex();

   if (fMCMCFlagsFillHistograms[index1] == false || fMCMCFlagsFillHistograms[index2] == false) {
// don't print any error message, should be done upstream
//      BCLog::OutError(
//            Form("BCModel::GetMarginalized : Distribuion for '%s' and/or '%s' not filled.",
//                  par1->GetName().data(), par2->GetName().data()));
      return 0;
   }

   if (index1 == index2) {
// don't print any error message, should be done upstream
//      BCLog::OutError("BCModel::GetMarginalized : Provided parameters are identical. Distribution not available.");
      return 0;
   }

   BCParameter * npar1 = par1;
   BCParameter * npar2 = par2;

   if (index1 > index2) {
      npar1 = par2;
      npar2 = par1;

      int itmp = index1;
      index1 = index2;
      index2 = itmp;
   }

   // get histogram
   TH2D * hist = MCMCGetH2Marginalized(index1, index2);

   if (hist == 0)
      return 0;

   BCH2D * hprob = new BCH2D();

   // set axis labels
   hist->SetName(Form("hist_%s_%s_%s", GetName().data(), npar1->GetName().data(), npar2->GetName().data()));
   hist->SetXTitle(Form("%s", npar1->GetName().data()));
   hist->SetYTitle(Form("%s", npar2->GetName().data()));
   hist->SetStats(kFALSE);

   double gmode[] = { fBestFitParameters.at(index1), fBestFitParameters.at(index2) };
   hprob->SetGlobalMode(gmode);

   // set histogram
   hprob->SetHistogram(hist);

   return hprob;
}

BCH2D* BCModel::GetMarginalized	(	const char *	name1,
		const char *	name2
	)			`[inline]`

Definition at line 614 of file BCModel.h.

         { return this -> GetMarginalized(this -> GetParameter(name1), this -> GetParameter(name2)); };

double BCModel::GetModelAPosterioriProbability ( ) [inline]

Returns:: The a posteriori probability.

Definition at line 96 of file BCModel.h.

         { return fModelAPosteriori; };

double BCModel::GetModelAPrioriProbability ( ) [inline]

Returns:: The a priori probability.

Definition at line 91 of file BCModel.h.

         { return fModelAPriori; };

std::string BCModel::GetName ( ) [inline]

Returns:: The name of the model.

Definition at line 81 of file BCModel.h.

         { return fName; };

int BCModel::GetNDataPoints ( )

Returns:: The number of data points in the current data set.

Definition at line 120 of file BCModel.cxx.

{
   int npoints = 0;
   if (fDataSet)
      npoints = fDataSet->GetNDataPoints();
   else {
      BCLog::OutWarning("BCModel::GetNDataPoints() : No data set defined.");
      return ERROR_NOEVENTS;
   }

   return npoints;
}

unsigned int BCModel::GetNDataPointsMaximum ( ) [inline]

Returns:: The maximum number of data points.

Definition at line 152 of file BCModel.h.

         { return fNDataPointsMaximum; };

unsigned int BCModel::GetNDataPointsMinimum ( ) [inline]

Returns:: The minimum number of data points.

Definition at line 147 of file BCModel.h.

         { return fNDataPointsMinimum; };

double BCModel::GetNormalization ( ) [inline]

Returns:: The normalization factor of the probability

Definition at line 101 of file BCModel.h.

         { return fNormalization; };

unsigned int BCModel::GetNParameters ( ) [inline]

Returns:: The number of parameters of the model.

Definition at line 157 of file BCModel.h.

         { return fParameterSet ? fParameterSet -> size() : 0; };

BCParameter * BCModel::GetParameter ( int index )

Parameters:

index

The index of the parameter in the parameter set.

Returns:: The parameter.

Definition at line 144 of file BCModel.cxx.

{
   if (!fParameterSet)
      return 0;

   if (index < 0 || index >= (int)GetNParameters()) {
      BCLog::OutWarning(
            Form("BCModel::GetParameter : Parameter index %d not within range.",index));
      return 0;
   }

   return fParameterSet->at(index);
}

BCParameter * BCModel::GetParameter ( const char * name )

Parameters:

name

The name of the parameter in the parameter set.

Returns:: The parameter.

Definition at line 159 of file BCModel.cxx.

{
   if (!fParameterSet)
      return 0;

   int index = -1;
   for (unsigned int i = 0; i < GetNParameters(); i++)
      if (name == GetParameter(i)->GetName())
         index = i;

   if (index < 0) {
      BCLog::OutWarning(
            Form("BCModel::GetParameter : Model %s has no parameter named '%s'",
                 GetName().data(), name));
      return 0;
   }

   return GetParameter(index);
}

BCParameterSet* BCModel::GetParameterSet ( ) [inline]

Returns:: parameter set

Definition at line 172 of file BCModel.h.

         { return fParameterSet; };

double BCModel::GetPValue ( ) [inline]

Definition at line 681 of file BCModel.h.

         { return fPValue; };

double BCModel::GetPvalueFromChi2	(	std::vector< double >	par,
		int	sigma_index
	)

Calculate p-value from Chi2 distribution for Gaussian problems

Parameters:

	par	Parameter set for the calculation of the likelihood
	sigma_index	Index of the sigma/uncertainty for the data points (for data in format "x y erry" the index would be 2)

Definition at line 1294 of file BCModel.cxx.

{
   double ll = LogLikelihood(par);
   int n = GetNDataPoints();

   double sum_sigma = 0;
   for (int i = 0; i < n; i++)
      sum_sigma += log(GetDataPoint(i)->GetValue(sigma_index));

   double chi2 = -2. * (ll + (double) n / 2. * log(2. * M_PI) + sum_sigma);

   fPValue = TMath::Prob(chi2, n);

   return fPValue;
}

double BCModel::GetPvalueFromChi2Johnson ( std::vector< double > par )

Calculate p-value from asymptotic Chi2 distribution for arbitrary problems using the definition (3) from Johnson, V.E. A Bayesian chi2 Test for Goodness-of-Fit. The Annals of Statistics 32, 2361-2384(2004).

Parameters:

par

Parameter set for the calculation of the likelihood

Definition at line 1318 of file BCModel.cxx.

{
   double chi2 = GetChi2Johnson(par);
   // look up corresponding p value
   fPValue = TMath::Prob(chi2, NumberBins() - 1);
   return fPValue;
}

double BCModel::GetPvalueFromChi2NDoF	(	std::vector< double >	par,
		int	sigma_index
	)

Definition at line 1493 of file BCModel.cxx.

{
   double ll = LogLikelihood(par);
   int n = GetNDataPoints();
   int npar = GetNParameters();

   double sum_sigma = 0;
   for (int i = 0; i < n; i++)
      sum_sigma += log(GetDataPoint(i)->GetValue(sigma_index));

   double chi2 = -2. * (ll + (double) n / 2. * log(2. * M_PI) + sum_sigma);

   fChi2NDoF = chi2 / double(n - npar);
   fPValueNDoF = TMath::Prob(chi2, n - npar);

   return fPValueNDoF;
}

double BCModel::GetPvalueFromKolmogorov	(	const std::vector< double > &	par,
		int	index
	)

Calculate p-value from Kolmogorov-Smirnov test statistic for 1D - datasets.

Parameters:

	par	Parameter set for the calculation of the likelihood
	index	Index of the data point in the BCDataSet (for data in format "x y erry" the index would be 1)

Definition at line 1512 of file BCModel.cxx.

{
   if (flag_discrete) {
      BCLog::OutError(Form("BCModel::GetPvalueFromKolmogorov : "
         "test defined only for continuous distributions."));
      return -1.;
   }

   //calculate the ECDF from the 1D data
   std::vector<double> yData = fDataSet->GetDataComponents(index);
   TH1D * ECDF = BCMath::ECDF(yData);

   int N = GetNDataPoints();

   // calculated expected CDF for unique points only
   std::set<double> uniqueObservations;
   for (int i = 0; i < N; i++)
       uniqueObservations.insert(CDF(par, i, false));

   int nUnique = uniqueObservations.size();
   if (nUnique != ECDF->GetNbinsX() + 1) {
      BCLog::OutError(Form("BCModel::GetPvalueFromKolmogorov : "
         "Number of unique data points doesn't match (%d vs %d)", nUnique,
            ECDF->GetNbinsX() + 1));
      return -1.;
   }

   // find maximum distance
   double distMax = 0.;

   // current distance
   double dist = 0.;

   std::set<double>::const_iterator iter = uniqueObservations.begin();
   for (int iBin = 0; iBin < nUnique; ++iBin) {
      // distance between current points
      dist = TMath::Abs(*iter - ECDF->GetBinContent(iBin + 1));
      // update maximum if necessary
      distMax = TMath::Max(dist, distMax);

//      BCLog::OutDebug(Form("BCModel::GetPvalueFromKolmogorov : "
//         "expected vs empirical (%f vs %f)", *iter, ECDF->GetBinContent(iBin
//            + 1)));

      // advance to next entry in the set
      ++iter;
   }

   // correct for total #points, not unique #points.
   // would need sqrt(n1*n2/(n1+n2)) if we had two experimental datasets
   double z = distMax * sqrt(N);

   fPValue = TMath::KolmogorovProb(z);

//   BCLog::OutDebug(Form("BCModel::GetPvalueFromKolmogorov : "
//      "max distance vs corrected (%f vs %f)", distMax, z));

   // clean up
   delete ECDF;

   return fPValue;
}

double BCModel::GetPValueNDoF ( ) [inline]

Definition at line 684 of file BCModel.h.

         { return fPValueNDoF; };

double BCModel::HessianMatrixElement	(	BCParameter *	parameter1,
		BCParameter *	parameter2,
		std::vector< double >	point
	)

Calculates the matrix element of the Hessian matrix

Parameters:

	parameter1	The parameter for the first derivative
	parameter2	The parameter for the first derivative

Returns:: The matrix element of the Hessian matrix

Definition at line 1626 of file BCModel.cxx.

{
   // check number of entries in vector
   if (point.size() != GetNParameters()) {
      BCLog::OutError("BCModel::HessianMatrixElement : Invalid number of entries in the vector.");
      return -1;
   }

   // define steps
   double nsteps = 1e5;
   double dx1 = par1->GetRangeWidth() / nsteps;
   double dx2 = par2->GetRangeWidth() / nsteps;

   // define points at which to evaluate
   std::vector<double> xpp = point;
   std::vector<double> xpm = point;
   std::vector<double> xmp = point;
   std::vector<double> xmm = point;

   int idx1 = par1->GetIndex();
   int idx2 = par2->GetIndex();

   xpp[idx1] += dx1;
   xpp[idx2] += dx2;

   xpm[idx1] += dx1;
   xpm[idx2] -= dx2;

   xmp[idx1] -= dx1;
   xmp[idx2] += dx2;

   xmm[idx1] -= dx1;
   xmm[idx2] -= dx2;

   // calculate probability at these points
   double ppp = Likelihood(xpp);
   double ppm = Likelihood(xpm);
   double pmp = Likelihood(xmp);
   double pmm = Likelihood(xmm);

   // return derivative
   return (ppp + pmm - ppm - pmp) / (4. * dx1 * dx2);
}

double BCModel::Likelihood ( std::vector< double > parameter ) [inline]

Returns the likelihood

Parameters:

parameters

A set of parameter values

Returns:: The likelihood

Definition at line 459 of file BCModel.h.

         { return exp( this->LogLikelihood(parameter) ); };

double BCModel::LogAPrioriProbability ( std::vector< double > parameters ) [virtual]

Returns natural logarithm of the prior probability. Method needs to be overloaded by the user.

Parameters:

parameters

A set of parameter values

Returns:: The prior probability p(parameters)

See also:: GetPrior(std::vector <double> parameters)

Reimplemented in BCGoFTest, BCSummaryPriorModel, and BCRooInterface.

Definition at line 570 of file BCModel.cxx.

{
   double logprob = 0;

   // get number of parameters
   int npar = GetNParameters();

   // loop over all 1-d priors
   for (int i = 0; i < npar; ++i) {
      if (fPriorContainer[i]) {
         // check what type of object is stored
         TF1 * f = dynamic_cast<TF1*>(fPriorContainer[i]);
         TH1 * h = dynamic_cast<TH1*>(fPriorContainer[i]);

         if (f) // TF1
            logprob += log(f->Eval(parameters[i]));
         else if (h) { // TH1
            if(fPriorContainerInterpolate[i])
               logprob += log(h->Interpolate(parameters[i]));
            else
               logprob += log(h->GetBinContent(h->FindBin(parameters[i])));
         }
         else
            BCLog::OutError(Form(
               "BCModel::LogAPrioriProbability : Prior for parameter %s "
               "is defined but not recodnized.",
               GetParameter(i)->GetName().c_str())); // this should never happen
      }
      // use constant only if user has defined it
      else if (!fPriorContainerConstant[i]) {
         BCLog::OutWarning(Form(
               "BCModel::LogAPrioriProbability: Prior for parameter %s "
               "is undefined. Using constant prior to proceed.",
               GetParameter(i)->GetName().c_str()));
         logprob -= log(GetParameter(i)->GetRangeWidth());
      }
   }

   // add the contribution from constant priors in one step
   logprob += fPriorConstantValue;

   return logprob;
}

virtual double BCModel::LogConditionalProbabilityEntry	(	BCDataPoint *	,
		std::vector< double >
	)			`[inline, virtual]`

Returns a natural logarithm of conditional probability. Method needs to be overloaded by the user.

Parameters:

	datapoint	A data point
	parameters	A set of parameter values

Returns:: The conditional probability p(datapoint|parameters)

See also:: GetConditionalEntry(BCDataPoint* datapoint, std::vector <double> parameters)

Definition at line 513 of file BCModel.h.

         { flag_ConditionalProbabilityEntry = false; return 0.; };

double BCModel::LogEval ( std::vector< double > parameters ) [virtual]

Overloaded function to evaluate integral.

Reimplemented from BCIntegrate.

Definition at line 615 of file BCModel.cxx.

{
   return LogProbabilityNN(parameters);
}

double BCModel::LogLikelihood ( std::vector< double > parameter ) [virtual]

Calculates natural logarithm of the likelihood. Method needs to be overloaded by the user.

Parameters:

parameters

A set of parameter values

Returns:: Natural logarithm of the likelihood

Reimplemented in BCGoFTest, BCSummaryPriorModel, BCEfficiencyFitter, BCGraphFitter, BCHistogramFitter, BCRooInterface, and BCTemplateFitter.

Definition at line 556 of file BCModel.cxx.

{
   double logprob = 0.;

   // add log of conditional probabilities event-by-event
   for (unsigned int i = 0; i < fDataSet->GetNDataPoints(); i++) {
      BCDataPoint * datapoint = GetDataPoint(i);
      logprob += LogConditionalProbabilityEntry(datapoint, parameters);
   }

   return logprob;
}

double BCModel::LogProbability ( std::vector< double > parameter )

Returns natural logarithm of the a posteriori probability given a set of parameter values

Parameters:

parameters

A set of parameter values

Returns:: The a posteriori probability

Definition at line 544 of file BCModel.cxx.

{
   // check if normalized
   if (fNormalization < 0. || fNormalization == 0.) {
      BCLog::OutError("BCModel::LogProbability. Normalization not available or zero.");
      return 0.;
   }

   return LogProbabilityNN(parameters) - log(fNormalization);
}

double BCModel::LogProbabilityNN ( std::vector< double > parameter )

Returns the natural logarithm of likelihood times prior probability given a set of parameter values

Parameters:

parameters

A set of parameter values

Returns:: The likelihood times prior probability

Definition at line 529 of file BCModel.cxx.

{
   // add log of conditional probability
   double logprob = LogLikelihood(parameters);

   // add log of prior probability
   if(fPriorConstantAll)
      logprob += fPriorConstantValue;
   else
      logprob += LogAPrioriProbability(parameters);

   return logprob;
}

int BCModel::MarginalizeAll ( )

Marginalize all probabilities wrt. single parameters and all combinations of two parameters. The individual distributions can be retrieved using the GetMarginalized method.

Returns:: Total number of marginalized distributions

Definition at line 822 of file BCModel.cxx.

{
   if(fParameterSet->size()<1) {
      BCLog::OutError(Form("MarginalizeAll : No parameters defined in model \'%s\'. Aborting.",GetName().data()));
      return 0;
   }

   BCLog::OutSummary(Form("Running MCMC for model \'%s\'",GetName().data()));

   // prepare function fitting
   double dx = 0.;
   double dy = 0.;

   if (fFitFunctionIndexX >= 0) {
      dx = (fDataPointUpperBoundaries->GetValue(fFitFunctionIndexX)
            - fDataPointLowerBoundaries->GetValue(fFitFunctionIndexX))
            / (double) fErrorBandNbinsX;

      dy = (fDataPointUpperBoundaries->GetValue(fFitFunctionIndexY)
            - fDataPointLowerBoundaries->GetValue(fFitFunctionIndexY))
            / (double) fErrorBandNbinsY;

      fErrorBandXY
            = new TH2D(TString::Format("errorbandxy_%d", BCLog::GetHIndex()), "",
                  fErrorBandNbinsX,
                  fDataPointLowerBoundaries->GetValue(fFitFunctionIndexX) - .5 * dx,
                  fDataPointUpperBoundaries->GetValue(fFitFunctionIndexX) + .5 * dx,
                  fErrorBandNbinsY,
                  fDataPointLowerBoundaries->GetValue(fFitFunctionIndexY) - .5 * dy,
                  fDataPointUpperBoundaries->GetValue(fFitFunctionIndexY) + .5 * dy);
      fErrorBandXY->SetStats(kFALSE);

      for (int ix = 1; ix <= fErrorBandNbinsX; ++ix)
         for (int iy = 1; iy <= fErrorBandNbinsX; ++iy)
            fErrorBandXY->SetBinContent(ix, iy, 0.);
   }

   MCMCMetropolis();
   FindModeMCMC();

   //   PrintResults(Form("%s.txt", GetName().data()));

   return 1;
}

double BCModel::Normalize ( )

Integrates over the un-normalized probability and updates fNormalization.

Definition at line 636 of file BCModel.cxx.

{
   if(fParameterSet->size()<1) {
      BCLog::OutError(Form("Normalize : No parameters defined in model \'%s\'. Aborting.",GetName().data()));
      return -1.;
   }

   BCLog::OutSummary(Form("Model \'%s\': Normalizing probability",GetName().data()));

   unsigned int n = GetNvar();

   // initialize BCIntegrate if not done already
   if (n == 0) {
      SetParameters(fParameterSet);
      n = GetNvar();
   }

   // integrate and get best fit parameters
   // maybe we have to remove the mode finding from here in the future
   fNormalization = Integrate();

   BCLog::OutDetail(Form(" --> Normalization factor : %.6g", fNormalization));

   return fNormalization;
}

int BCModel::NumberBins ( ) [inline, private]

rule of thumb for good number of bins (Wald1942, Johnson2004) to group observations updated so minimum is three bins (for 1-5 observations)!

Parameters:

Definition at line 866 of file BCModel.h.

         { return (int)(exp(0.4 * log(this -> GetNDataPoints())) + 2); }

int BCModel::PrintAllMarginalized	(	const char *	file,
		unsigned int	hdiv = `1`,
		unsigned int	ndiv = `1`
	)

Definition at line 1053 of file BCModel.cxx.

{
   if (!fMCMCFlagFillHistograms) {
      BCLog::OutError("BCModel::PrintAllMarginalized : Marginalized distributions not filled.");
      return 0;
   }

   int npar = GetNParameters();

   if (fMCMCH1Marginalized.size() == 0 || (fMCMCH2Marginalized.size() == 0
         && npar > 1)) {
      BCLog::OutError(
            "BCModel::PrintAllMarginalized : Marginalized distributions not available.");
      return 0;
   }

   // if there's only one parameter, we just want to call Print()
   if (fMCMCH1Marginalized.size() == 1 && fMCMCH2Marginalized.size() == 0) {
      BCParameter * a = GetParameter(0);
      if (GetMarginalized(a))
         GetMarginalized(a)->Print(file);
      return 1;
   }

   int c_width = 600; // default canvas width
   int c_height = 600; // default canvas height

   int type = 112; // landscape

   if (hdiv > vdiv) {
      if (hdiv > 3) {
         c_width = 1000;
         c_height = 700;
      }
      else {
         c_width = 800;
         c_height = 600;
      }
   }
   else if (hdiv < vdiv) {
      if (hdiv > 3) {
         c_height = 1000;
         c_width = 700;
      }
      else {
         c_height = 800;
         c_width = 600;
      }
      type = 111;
   }

   // calculate number of plots
   int nplots2d = npar * (npar - 1) / 2;
   int nplots = npar + nplots2d;

   // give out warning if too many plots
   BCLog::OutSummary(
         Form(
               "Printing all marginalized distributions (%d x 1D + %d x 2D = %d) into file %s",
               npar, nplots2d, nplots, file));
   if (nplots > 100)
      BCLog::OutDetail("This can take a while...");

   // setup the canvas and postscript file
   TCanvas * c = new TCanvas("c", "canvas", c_width, c_height);

   TPostScript * ps = new TPostScript(file, type);

   if (type == 112)
      ps->Range(24, 16);
   else
      ps->Range(16, 24);

   // draw all 1D distributions
   ps->NewPage();
   c->cd();
   c->Clear();
   c->Divide(hdiv, vdiv);

   int n = 0;
   for (int i = 0; i < npar; i++) {
      // get corresponding parameter
      BCParameter * a = GetParameter(i);

      // check if histogram exists
      if (!GetMarginalized(a))
         continue;

      // check if histogram is filled
      if (GetMarginalized(a)->GetHistogram()->Integral() <= 0)
         continue;

      // if current page is full, switch to new page
      if (i != 0 && i % (hdiv * vdiv) == 0) {
         c->Update();
         ps->NewPage();
         c->cd();
         c->Clear();
         c->Divide(hdiv, vdiv);
      }

      // go to next pad
      c->cd(i % (hdiv * vdiv) + 1);

      GetMarginalized(a)->Draw();
      n++;

      if (n % 100 == 0)
         BCLog::OutDetail(Form(" --> %d plots done", n));
   }

   if (n > 0) {
      c->Update();

      // draw all the 2D distributions
      ps->NewPage();
      c->cd();
      c->Clear();
   }

   c->Divide(hdiv, vdiv);

   int k = 0;
   for (int i = 0; i < npar - 1; i++) {
      if (!fMCMCFlagsFillHistograms[i])
         continue;
      for (int j = i + 1; j < npar; j++) {
         if (!fMCMCFlagsFillHistograms[j])
            continue;

         // get corresponding parameters
         BCParameter * a = GetParameter(i);
         BCParameter * b = GetParameter(j);

         // check if histogram exists
         if (!GetMarginalized(a, b))
            continue;

         // check if histogram is filled
         if (GetMarginalized(a, b)->GetHistogram()->Integral() <= 0)
            continue;

         // if current page is full, switch to new page
         if (k != 0 && k % (hdiv * vdiv) == 0) {
            c->Update();
            ps->NewPage();
            c->cd();
            c->Clear();
            c->Divide(hdiv, vdiv);
         }

         // go to next pad
         c->cd(k % (hdiv * vdiv) + 1);

         double meana = (a->GetLowerLimit() + a->GetUpperLimit()) / 2.;
         double deltaa = (a->GetUpperLimit() - a->GetLowerLimit());
         if (deltaa <= 1e-7 * meana)
            continue;

         double meanb = (b->GetLowerLimit() + b->GetUpperLimit()) / 2.;
         double deltab = (b->GetUpperLimit() - b->GetLowerLimit());
         if (deltab <= 1e-7 * meanb)
            continue;

         GetMarginalized(a, b)->Draw(52);
         k++;

         if ((n + k) % 100 == 0)
            BCLog::OutDetail(Form(" --> %d plots done", n + k));
      }
   }

   if ((n + k) > 100 && (n + k) % 100 != 0)
      BCLog::OutDetail(Form(" --> %d plots done", n + k));

   c->Update();
   ps->Close();

   delete c;
   delete ps;

   // return total number of drawn histograms
   return n + k;
}

int BCModel::PrintAllMarginalized1D ( const char * filebase )

Definition at line 1000 of file BCModel.cxx.

{
   if (fMCMCH1Marginalized.size() == 0) {
      BCLog::OutError("BCModel::PrintAllMarginalized : Marginalized distributions not available.");
      return 0;
   }

   int n = GetNParameters();
   for (int i = 0; i < n; i++) {
      BCParameter * a = GetParameter(i);
      if (GetMarginalized(a))
         GetMarginalized(a)->Print(Form("%s_1D_%s.ps", filebase, a->GetName().data()));
   }

   return n;
}

int BCModel::PrintAllMarginalized2D ( const char * filebase )

Definition at line 1018 of file BCModel.cxx.

{
   if (fMCMCH2Marginalized.size() == 0) {
      BCLog::OutError("BCModel::PrintAllMarginalized : Marginalized distributions not available.");
      return 0;
   }

   int k = 0;
   int n = GetNParameters();
   for (int i = 0; i < n - 1; i++) {
      for (int j = i + 1; j < n; j++) {
         BCParameter * a = GetParameter(i);
         BCParameter * b = GetParameter(j);

         double meana = (a->GetLowerLimit() + a->GetUpperLimit()) / 2.;
         double deltaa = (a->GetUpperLimit() - a->GetLowerLimit());
         if (deltaa <= 1e-7 * meana)
            continue;

         double meanb = (b->GetLowerLimit() + b->GetUpperLimit()) / 2.;
         double deltab = (b->GetUpperLimit() - b->GetLowerLimit());
         if (deltab <= 1e-7 * meanb)
            continue;

         if (GetMarginalized(a, b))
            GetMarginalized(a, b)->Print(
                  Form("%s_2D_%s_%s.ps",filebase, a->GetName().data(), b->GetName().data()));
         k++;
      }
   }

   return k;
}

void BCModel::PrintHessianMatrix ( std::vector< double > parameters )

Prints matrix elements of the Hessian matrix

Parameters:

parameters

The parameter values at which point to evaluate the matrix

Definition at line 2249 of file BCModel.cxx.

{
   // check number of entries in vector
   if (parameters.size() != GetNParameters()) {
      BCLog::OutError("BCModel::PrintHessianMatrix : Invalid number of entries in the vector");
      return;
   }

   // print to screen
   std::cout << std::endl << " Hessian matrix elements: " << std::endl << " Point: ";

   for (int i = 0; i < int(parameters.size()); i++)
      std::cout << parameters.at(i) << " ";
   std::cout << std::endl;

   // loop over all parameter pairs
   for (unsigned int i = 0; i < GetNParameters(); i++)
      for (unsigned int j = 0; j < i; j++) {
         // calculate Hessian matrix element
         double hessianmatrixelement = HessianMatrixElement(
               fParameterSet->at(i), fParameterSet->at(j), parameters);

         // print to screen
         std::cout << " " << i << " " << j << " : " << hessianmatrixelement << std::endl;
      }
}

void BCModel::PrintResults ( const char * file )

Prints a summary of the Markov Chain Monte Carlo to a file.

Definition at line 2042 of file BCModel.cxx.

{
   // print summary of Markov Chain Monte Carlo

   // open file
   ofstream ofi(file);

   // check if file is open
   if (!ofi.is_open()) {
      std::cerr << "Couldn't open file " << file << std::endl;
      return;
   }

   // number of parameters and chains
   int npar = MCMCGetNParameters();
   int nchains = MCMCGetNChains();

   // check convergence
   bool flag_conv = MCMCGetNIterationsConvergenceGlobal() > 0;

   ofi << std::endl 
          << " -----------------------------------------------------" << std::endl 
          << " Summary" << std::endl
       << " -----------------------------------------------------" << std::endl
       << std::endl;

   ofi << " Model summary" << std::endl << " =============" << std::endl
       << " Model: " << fName.data() << std::endl
       << " Number of parameters: " << npar << std::endl
       << " List of Parameters and ranges:" << std::endl;
   for (int i = 0; i < npar; ++i)
      ofi << "  (" << i << ") Parameter \""
          << fParameterSet->at(i)->GetName().data() << "\"" << ": "
          << "(" << fParameterSet->at(i)->GetLowerLimit() << ", "
          << fParameterSet->at(i)->GetUpperLimit() << ")" << std::endl;
   ofi << std::endl;

   // give warning if MCMC did not converge
   if (!flag_conv && fMCMCFlagRun)
      ofi << " WARNING: the Markov Chain did not converge!" << std::endl
          << " Be cautious using the following results!" << std::endl
          << std::endl;

   // print results of marginalization (if MCMC was run)
   if (fMCMCFlagRun && fMCMCFlagFillHistograms) {
      ofi << " Results of the marginalization" << std::endl
          << " ==============================" << std::endl
          << " List of parameters and properties of the marginalized"
          << std::endl << " distributions:" << std::endl;
      for (int i = 0; i < npar; ++i) {
         if (!fMCMCFlagsFillHistograms[i])
            continue;

         BCH1D * bch1d = GetMarginalized(fParameterSet->at(i));

         ofi << "  (" << i << ") Parameter \""
             << fParameterSet->at(i)->GetName().data() << "\":" << std::endl

             << "      Mean +- sqrt(V):                " << std::setprecision(4)
             << bch1d->GetMean() << " +- " << std::setprecision(4)
             << bch1d->GetRMS() << std::endl 

                   << "      Median +- central 68% interval: "
             << std::setprecision(4) << bch1d->GetMedian() << " +  "
             << std::setprecision(4) << bch1d->GetQuantile(0.84) - bch1d->GetMedian()
             << " - " << std::setprecision(4)
             << bch1d->GetMedian() - bch1d->GetQuantile(0.16) << std::endl

             << "      (Marginalized) mode:            " << bch1d->GetMode() << std::endl;

         ofi << "       5% quantile:                   " << std::setprecision(4)
             << bch1d->GetQuantile(0.05) << std::endl
             << "      10% quantile:                   " << std::setprecision(4)
             << bch1d->GetQuantile(0.10) << std::endl
             << "      16% quantile:                   " << std::setprecision(4)
             << bch1d->GetQuantile(0.16) << std::endl
             << "      84% quantile:                   " << std::setprecision(4)
             << bch1d->GetQuantile(0.85) << std::endl
             << "      90% quantile:                   " << std::setprecision(4)
             << bch1d->GetQuantile(0.90) << std::endl
             << "      95% quantile:                   " << std::setprecision(4)
             << bch1d->GetQuantile(0.95) << std::endl;

         std::vector<double> v;
         v = bch1d->GetSmallestIntervals(0.68);
         int ninter = int(v.size());
             
             ofi << "      Smallest interval(s) containing 68% and local modes:"
                   << std::endl;
         for (int j = 0; j < ninter; j += 5)
                ofi << "       (" << v[j] << ", " << v[j + 1]
                      << ") (local mode at " << v[j + 3] << " with rel. height "
                      << v[j + 2] << "; rel. area " << v[j + 4] << ")"
                      << std::endl;
             ofi << std::endl;
      }
   }

   ofi << " Results of the optimization" << std::endl
       << " ===========================" << std::endl
       << " Optimization algorithm used:"
       << DumpOptimizationMethod() << std::endl;

   if (int(fBestFitParameters.size()) > 0) {
      ofi << " List of parameters and global mode:" << std::endl;
      for (int i = 0; i < npar; ++i) {
         ofi << "  (" << i << ") Parameter \""
             << fParameterSet->at(i)->GetName().data() << "\": "
             << fBestFitParameters[i];
             if (int(fBestFitParameterErrors.size()) == npar)
                if(fBestFitParameterErrors[i]>=0.)
                   ofi << " +- " << fBestFitParameterErrors[i];
         ofi << std::endl;
      }
      ofi << std::endl;
   }
   else {
      ofi << " No best fit information available." << std::endl;
      ofi << std::endl;
   }

   if (fPValue >= 0.) {
      ofi << " Results of the model test" << std::endl
          << " =========================" << std::endl
          << " p-value at global mode: " << fPValue << std::endl << std::endl;
   }

   if (fNormalization >= 0.) {
       ofi << " Results of the normalization" << std::endl
             << " ============================" << std::endl
             << " Integration method used:"
             << DumpIntegrationMethod() << std::endl;
       ofi << " Normalization factor: " << fNormalization << std::endl << std::endl;
   }

   if (fMCMCFlagRun) {
      ofi << " Status of the MCMC" << std::endl << " ==================" << std::endl
               << " Convergence reached:                    " << (flag_conv ? "yes" : "no")
          << std::endl;

      if (flag_conv)
         ofi << " Number of iterations until convergence: "
             << MCMCGetNIterationsConvergenceGlobal() << std::endl;
      else
         ofi << " WARNING: the Markov Chain did not converge! Be\n"
             << " cautious using the following results!" << std::endl
             << std::endl;
      ofi << " Number of chains:                       " << MCMCGetNChains() << std::endl
          << " Number of iterations per chain:         " << MCMCGetNIterationsRun() << std::endl
          << " Average efficiencies:" << std::endl;

      std::vector<double> efficiencies;
      efficiencies.assign(npar, 0.);

      for (int ipar = 0; ipar < npar; ++ipar)
         for (int ichain = 0; ichain < nchains; ++ichain) {
            int index = ichain * npar + ipar;
            efficiencies[ipar] +=
                  double(MCMCGetNTrialsTrue().at(index)) / double(MCMCGetNTrialsTrue().at(index)
                  + MCMCGetNTrialsFalse().at(index)) / double(nchains) * 100.;
         }

      for (int ipar = 0; ipar < npar; ++ipar)
         ofi << "  (" << ipar << ") Parameter \""
             << fParameterSet->at(ipar)->GetName().data() << "\": "
             << efficiencies.at(ipar) << "%" << std::endl;
      ofi << std::endl;
   }

    ofi << " -----------------------------------------------------" << std::endl
          << " Notation: " << std::endl
          << " Mean        : mean value of the marg. pdf" << std::endl
          << " Median      : maximum of the marg. pdf" << std::endl
          << " Marg. mode  : most probable value of the marg. pdf" << std::endl
          << " V           : Variance of the marg. pdf" << std::endl
          << " Quantiles   : most commonly used quantiles" <<std::endl
          << " -----------------------------------------------------" << std::endl
          << std::endl;

   // close file
   //   ofi.close;
}

void BCModel::PrintShortFitSummary ( int chi2flag = 0 )

Prints a short summary of the fit results on the screen.

Definition at line 2226 of file BCModel.cxx.

{
   BCLog::OutSummary("---------------------------------------------------");
   BCLog::OutSummary(Form("Fit summary for model \'%s\':", GetName().data()));
   BCLog::OutSummary(Form("   Number of parameters:  Npar  = %i", GetNParameters()));
   BCLog::OutSummary(Form("   Number of data points: Ndata = %i", GetNDataPoints()));
   BCLog::OutSummary("   Number of degrees of freedom:");
   BCLog::OutSummary(Form("      NDoF = Ndata - Npar = %i", GetNDataPoints() - GetNParameters()));

   BCLog::OutSummary("   Best fit parameters (global):");
   for (unsigned int i = 0; i < GetNParameters(); ++i)
      BCLog::OutSummary(Form("      %s = %.3g", GetParameter(i)->GetName().data(), GetBestFitParameter(i)));

   BCLog::OutSummary("   Goodness-of-fit test:");
   BCLog::OutSummary(Form("      p-value = %.3g", GetPValue()));
   if (chi2flag) {
      BCLog::OutSummary(Form("      p-value corrected for NDoF = %.3g", GetPValueNDoF()));
      BCLog::OutSummary(Form("      chi2 / NDoF = %.3g", GetChi2NDoF()));
   }
   BCLog::OutSummary("---------------------------------------------------");
}

void BCModel::PrintSummary ( )

Prints a summary on the screen.

Definition at line 1994 of file BCModel.cxx.

{
   int nparameters = GetNParameters();

   // model summary
   std::cout << std::endl
             << "   ---------------------------------" << std::endl
             << "    Model : " << fName.data() << std::endl
             << "   ---------------------------------" << std::endl
             << "     Index                : " << fIndex << std::endl
             << "     Number of parameters : " << nparameters << std::endl << std::endl
             << "     - Parameters : " << std::endl << std::endl;

   // parameter summary
   for (int i = 0; i < nparameters; i++)
      fParameterSet->at(i)->PrintSummary();

   // best fit parameters
   if (GetBestFitParameters().size() > 0) {
      std::cout << std::endl << "     - Best fit parameters :" << std::endl << std::endl;

      for (int i = 0; i < nparameters; i++) {
         std::cout << "       " << fParameterSet->at(i)->GetName().data()
                   << " = " << GetBestFitParameter(i) << " (overall)" << std::endl;
         if ((int) fBestFitParametersMarginalized.size() == nparameters)
            std::cout << "       "
                      << fParameterSet->at(i)->GetName().data() << " = "
                      << GetBestFitParameterMarginalized(i)
                      << " (marginalized)" << std::endl;
      }
   }

   std::cout << std::endl;

   // model testing
   if (fPValue >= 0) {
      double likelihood = Likelihood(GetBestFitParameters());
      std::cout << "   - Model testing:" << std::endl << std::endl
                << "       p(data|lambda*) = " << likelihood << std::endl
                << "       p-value         = " << fPValue << std::endl << std::endl;
   }

   // normalization
   if (fNormalization > 0)
      std::cout << "     Normalization        : " << fNormalization << std::endl;
}

double BCModel::Probability ( std::vector< double > parameter ) [inline]

Returns the a posteriori probability given a set of parameter values

Parameters:

parameters

A set of parameter values

Returns:: The a posteriori probability

Definition at line 487 of file BCModel.h.

         { return exp( this->LogProbability(parameter) ); };

double BCModel::ProbabilityNN ( std::vector< double > parameter ) [inline]

Returns the likelihood times prior probability given a set of parameter values

Parameters:

parameters

A set of parameter values

Returns:: The likelihood times prior probability

Definition at line 473 of file BCModel.h.

         { return exp( this->LogProbabilityNN(parameter) ); };

int BCModel::ReadErrorBandFromFile ( const char * file )

Read

Definition at line 973 of file BCModel.cxx.

{
   TFile * froot = new TFile(file);
   if (!froot->IsOpen()) {
      BCLog::OutError(Form("BCModel::ReadErrorBandFromFile. Couldn't open file %s.", file));
      return 0;
   }

   int r = 0;

   TH2D * h2 = (TH2D*) froot->Get("errorbandxy");
   if (h2) {
      h2->SetDirectory(0);
      h2->SetName(TString::Format("errorbandxy_%d", BCLog::GetHIndex()));
      SetErrorBandHisto(h2);
      r = 1;
   }
   else
      BCLog::OutWarning(
            Form("BCModel::ReadErrorBandFromFile : Couldn't read histogram \"errorbandxy\" from file %s.",file));

   froot->Close();

   return r;
}

int BCModel::ReadMarginalizedFromFile ( const char * file )

Read

Definition at line 915 of file BCModel.cxx.

{
   TFile * froot = new TFile(file);
   if (!froot->IsOpen()) {
      BCLog::OutError(Form("BCModel::ReadMarginalizedFromFile : Couldn't open file %s.", file));
      return 0;
   }

   // We reset the MCMCEngine here for the moment.
   // In the future maybe we only want to do this if the engine
   // wans't initialized at all or when there were some changes
   // in the model.
   // But maybe we want reset everything since we're overwriting
   // the marginalized distributions anyway.
   MCMCInitialize();

   int k = 0;
   int n = GetNParameters();
   for (int i = 0; i < n; i++) {
      BCParameter * a = GetParameter(i);
      TKey * key = froot->GetKey(TString::Format("hist_%s_%s",GetName().data(), a->GetName().data()));
      if (key) {
         TH1D * h1 = (TH1D*) key->ReadObjectAny(TH1D::Class());
         h1->SetDirectory(0);
         if (SetMarginalized(i, h1))
            k++;
      }
      else
         BCLog::OutWarning(
               Form("BCModel::ReadMarginalizedFromFile : Couldn't read histogram \"hist_%s_%s\" from file %s.",
                     GetName().data(), a->GetName().data(), file));
   }

   for (int i = 0; i < n - 1; i++) {
      for (int j = i + 1; j < n; j++) {
         BCParameter * a = GetParameter(i);
         BCParameter * b = GetParameter(j);
         TKey * key = froot->GetKey(
               TString::Format("hist_%s_%s_%s",GetName().data(), a->GetName().data(),b->GetName().data()));
         if (key) {
            TH2D * h2 = (TH2D*) key->ReadObjectAny(TH2D::Class());
            h2->SetDirectory(0);
            if (SetMarginalized(i, j, h2))
               k++;
         }
         else
            BCLog::OutWarning(
                  Form("BCModel::ReadMarginalizedFromFile : Couldn't read histogram \"hist_%s_%s_%s\" from file %s.",
                        GetName().data(), a->GetName().data(), b->GetName().data(), file));
      }
   }

   froot->Close();

   return k;
}

int BCModel::ReadMode ( const char * file )

Read mode from file created by WriteMode() call

Definition at line 783 of file BCModel.cxx.

{
   ifstream ifi(file);
   if (!ifi.is_open()) {
      BCLog::OutError(Form("BCModel::ReadMode : Couldn't open file %s.", file));
      return 0;
   }

   int npar = fParameterSet->size();
   std::vector<double> mode;

   int i = 0;
   while (i < npar && !ifi.eof()) {
      double a;
      ifi >> a;
      mode.push_back(a);
      i++;
   }

   if (i < npar) {
      BCLog::OutError(Form("BCModel::ReadMode : Couldn't read mode from file %s.", file));
      BCLog::OutError(Form("BCModel::ReadMode : Expected %d parameters, found %d.", npar, i));
      return 0;
   }

   BCLog::OutSummary(
         Form("#  Read in best fit parameters (mode) for model \'%s\' from file %s:",
               fName.data(), file));
   SetMode(mode);
   for (int j = 0; j < npar; j++)
      BCLog::OutSummary(Form("#   ->Parameter %d : %s = %e", j,
            fParameterSet->at(j)->GetName().data(), fBestFitParameters[j]));

   BCLog::OutWarning("#  ! Best fit values obtained before this call will be overwritten !");

   return npar;
}

void BCModel::RecalculatePriorConstant ( ) [private]

Update the constant part of the prior.

Definition at line 1939 of file BCModel.cxx.

{
   fPriorConstantValue = 0.;

   // get number of parameters
   int npar = GetNParameters();

   int nconstant = 0;

   for (int i=0; i<npar; ++i)
      if (fPriorContainerConstant[i]) {
         fPriorConstantValue -= log(GetParameter(i)->GetRangeWidth());
         ++nconstant;
      }

   if (nconstant == npar)
      fPriorConstantAll = true;
   else
      fPriorConstantAll = false;
}

int BCModel::ResetResults ( )

Reset all results.

Returns:: An error code.

Definition at line 1983 of file BCModel.cxx.

{
   BCIntegrate::IntegrateResetResults();
  
   BCEngineMCMC::MCMCResetResults();

   // no error
   return 1;
}

double BCModel::SamplingFunction ( std::vector< double > parameters ) [virtual]

Sampling function used for importance sampling. Method needs to be overloaded by the user.

Parameters:

parameters

A set of parameter values

Returns:: The probability density at the parameter values

Definition at line 627 of file BCModel.cxx.

{
   double probability = 1.;
   for (std::vector<BCParameter *>::const_iterator it = fParameterSet->begin(); it != fParameterSet->end(); ++it)
      probability *= 1. / ((*it)->GetUpperLimit() - (*it)->GetLowerLimit());
   return probability;
}

void BCModel::SetDataBoundaries	(	unsigned int	index,
		double	lowerboundary,
		double	upperboundary,
		bool	fixed = `false`
	)

Definition at line 422 of file BCModel.cxx.

{
   // check if data set exists
   if (!fDataSet) {
      BCLog::OutError("BCModel::SetDataBoundaries : Need to define data set first.");
      return;
   }

   // check if index is within range
   if (index < 0 || index > fDataSet->GetDataPoint(0)->GetNValues()) {
      BCLog::OutError("BCModel::SetDataBoundaries : Index out of range.");
      return;
   }

   // check if boundary data points exist
   if (!fDataPointLowerBoundaries)
      fDataPointLowerBoundaries = new BCDataPoint(fDataSet->GetDataPoint(0)->GetNValues());

   if (!fDataPointUpperBoundaries)
      fDataPointUpperBoundaries = new BCDataPoint(fDataSet->GetDataPoint(0)->GetNValues());

   if (fDataFixedValues.size() == 0)
      fDataFixedValues.assign(fDataSet->GetDataPoint(0)->GetNValues(), false);

   // set boundaries
   fDataPointLowerBoundaries->SetValue(index, lowerboundary);
   fDataPointUpperBoundaries->SetValue(index, upperboundary);
   fDataFixedValues[index] = fixed;
}

void BCModel::SetDataSet ( BCDataSet * dataset ) [inline]

Sets the data set.

Parameters:

dataset

A data set

Definition at line 292 of file BCModel.h.

         { fDataSet = dataset; fNormalization = -1.0; };

void BCModel::SetErrorBandContinuous ( bool flag )

Sets the error band flag to continuous function

Definition at line 453 of file BCModel.cxx.

{
   fErrorBandContinuous = flag;

   if (flag)
      return;

   // clear x-values
   fErrorBandX.clear();

   // copy data x-values
   for (unsigned int i = 0; i < fDataSet->GetNDataPoints(); ++i)
      fErrorBandX.push_back(fDataSet->GetDataPoint(i)->GetValue(fFitFunctionIndexX));
}

void BCModel::SetGoFNChains ( int n ) [inline]

Definition at line 712 of file BCModel.h.

         { fGoFNChains=n; };

void BCModel::SetGoFNIterationsMax ( int n ) [inline]

Definition at line 700 of file BCModel.h.

         { fGoFNIterationsMax=n; };

void BCModel::SetGoFNIterationsRun ( int n ) [inline]

Definition at line 706 of file BCModel.h.

         { fGoFNIterationsRun=n; };

void BCModel::SetIndex ( int index ) [inline]

Sets the index of the model within the BCModelManager.

Parameters:

index

The index of the model

Definition at line 251 of file BCModel.h.

         { fIndex = index; };

void BCModel::SetModelAPosterioriProbability ( double probability ) [inline]

Sets the a posteriori probability for a model.

Parameters:

	model	The model
	probability	The a posteriori probability

Definition at line 279 of file BCModel.h.

         { fModelAPosteriori = probability; };

void BCModel::SetModelAPrioriProbability ( double probability ) [inline]

Sets the a priori probability for a model.

Parameters:

	model	The model
	probability	The a priori probability

Definition at line 272 of file BCModel.h.

         { fModelAPriori = probability; };

void BCModel::SetName ( const char * name ) [inline]

Sets the name of the model.

Parameters:

name

Name of the model

Definition at line 245 of file BCModel.h.

         { fName = name; };

void BCModel::SetNbins	(	const char *	parname,
		int	nbins
	)

Set the number of bins for the marginalized distribution of a parameter.

Parameters:

	parname	The name of the parameter in the parameter set
	nbins	Number of bins (default = 100)

Definition at line 231 of file BCModel.cxx.

{
   BCParameter * p = GetParameter(parname);
   if (!p) {
      BCLog::OutWarning(
            Form("BCModel::SetNbins : Parameter '%s' not found so Nbins not set",parname));
      return;
   }

   BCIntegrate::SetNbins(nbins, p->GetIndex());
}

void BCModel::SetNDataPointsMaximum ( unsigned int maximum ) [inline]

Sets the maximum number of data points.

Definition at line 309 of file BCModel.h.

         { fNDataPointsMaximum = maximum; };

void BCModel::SetNDataPointsMinimum ( unsigned int minimum ) [inline]

Sets the minimum number of data points.

Definition at line 304 of file BCModel.h.

         { fNDataPointsMinimum = minimum; };

void BCModel::SetNormalization ( double norm ) [inline]

Sets the normalization of the likelihood. The normalization is the integral of likelihood over all parameters.

Parameters:

norm

The normalization of the likelihood

Definition at line 286 of file BCModel.h.

         { fNormalization = norm; };

int BCModel::SetParameterRange	(	int	index,
		double	parmin,
		double	parmax
	)

Set the range of a parameter

Parameters:

	index	The parameter index
	parmin	The parameter minimum
	parmax	The parameter maximum

Returns:: An error code.

Definition at line 1961 of file BCModel.cxx.

{
   // check index
   if (index < 0 || index >= int(GetNParameters())) {
      BCLog::OutError("BCModel::SetParameterRange : Index out of range.");
      return 0; 
   }

   // set parameter ranges in BAT 
   GetParameter(index)->SetLowerLimit(parmin); 
   GetParameter(index)->SetUpperLimit(parmax); 
   fMCMCBoundaryMin[index] = parmin; 
   fMCMCBoundaryMax[index] = parmax; 

   // reset results
   ResetResults();

   // no error 
   return 1; 
}

void BCModel::SetParameterSet ( BCParameterSet * parset ) [inline]

Set all parameters of the model using a BCParameterSet container.

pointer to parameter set

Definition at line 257 of file BCModel.h.

         { fParameterSet = parset; };

int BCModel::SetPrior	(	const char *	name,
		TH1 *	h,
		bool	flag = `false`
	)

Set prior for a parameter.

Parameters:

	name	parameter name
	h	pointer to a histogram describing the prior
	flag	whether or not to use linear interpolation

Returns:: An error code.

Definition at line 1852 of file BCModel.cxx.

{
   // find index
   int index = -1;
   for (unsigned int i = 0; i < GetNParameters(); i++)
      if (name == GetParameter(i)->GetName())
        index = i;

   // set prior
   return SetPrior(index, h, interpolate);
}

int BCModel::SetPrior	(	const char *	name,
		TF1 *	f
	)

Set prior for a parameter.

Parameters:

	name	The parameter name
	f	A pointer to a function describing the prior

Returns:: An error code.

Definition at line 1725 of file BCModel.cxx.

{
   // find index
   int index = -1;
   for (unsigned int i = 0; i < GetNParameters(); i++)
      if (name == GetParameter(i)->GetName())
        index = i;

   // set prior
   return SetPrior(index, f);
}

int BCModel::SetPrior	(	int	index,
		TH1 *	h,
		bool	flag = `false`
	)

Set prior for a parameter.

Parameters:

	index	parameter index
	h	pointer to a histogram describing the prior
	flag	whether or not to use linear interpolation

Returns:: An error code.

Definition at line 1810 of file BCModel.cxx.

{
   // check index range
   if (index < 0 && index >= int(GetNParameters())) {
      BCLog::OutError("BCModel::SetPrior : Index out of range."); 
      return 0;
   }

   // if the histogram exists
   if(h) {

      // check if histogram is 1d
      if (h->GetDimension() != 1) {
         BCLog::OutError(Form("BCModel::SetPrior : Histogram given for parameter %d is not 1D.",index)); 
         return 0;
      }

      // normalize the histogram
      h->Scale(1./h->Integral("width"));

      if(fPriorContainer[index])
         delete fPriorContainer[index];

      // set function
      fPriorContainer[index] = new TH1(*h);

      if (interpolate)
         fPriorContainerInterpolate[index] = true;

      fPriorContainerConstant[index] = false;
   }

   RecalculatePriorConstant();

   // reset all results
   ResetResults();

   // no error 
   return 1; 
}

int BCModel::SetPrior	(	int	index,
		TF1 *	f
	)

Set prior for a parameter.

Parameters:

	index	The parameter index
	f	A pointer to a function describing the prior

Returns:: An error code.

Definition at line 1699 of file BCModel.cxx.

{
   // check index range
   if (index < 0 && index >= int(GetNParameters())) {
      BCLog::OutError("BCModel::SetPrior : Index out of range."); 
      return 0;
   }

   if (fPriorContainer[index])
      delete fPriorContainer[index];

   // copy function
   fPriorContainer[index] = new TF1(*f);

   fPriorContainerConstant[index] = false;

   RecalculatePriorConstant();

   // reset all results
   ResetResults();

   // no error 
   return 1; 
}

int BCModel::SetPriorConstant ( const char * name )

Set constant prior for this parameter

Parameters:

name

the name of the parameter

Returns:: An error code

Definition at line 1891 of file BCModel.cxx.

{
   // find index
   int index = -1;
   for (unsigned int i = 0; i < GetNParameters(); i++) {
      if (name == GetParameter(i)->GetName()) {
         index = i;
         break;
      }
   }

   if (index == -1) {
      BCLog::OutError(Form(
            "BCModel::SetPriorConstant : parameter '%s' doesn't exist.", name));
      return 0;
   }

   return SetPriorConstant(index);
}

int BCModel::SetPriorConstant ( int index )

Set constant prior for this parameter

Parameters:

index

the index of the parameter

Returns:: An error code

Definition at line 1865 of file BCModel.cxx.

{
   // check index range
   if (index < 0 && index >= int(GetNParameters())) {
      BCLog::OutError("BCModel::SetPriorConstant : Index out of range.");
      return 0;
   }

   if(fPriorContainer[index]) {
      delete fPriorContainer[index];
      fPriorContainer[index] = 0;
   }

   // set prior to a constant
   fPriorContainerConstant[index] = true;

   RecalculatePriorConstant();

   // reset all results
   ResetResults();

   // no error
   return 1;
}

int BCModel::SetPriorConstantAll ( )

Enable caching the constant value of the prior, so LogAPrioriProbability is called only once. Note that the prior for ALL parameters is assumed to be constant. The value is computed from the parameter ranges, so make sure these are defined before this method is called.

Returns:: An error code

Definition at line 1912 of file BCModel.cxx.

{
   // get number of parameters
   int nPar = GetNParameters();

   if (nPar == 0)
      BCLog::OutWarning("BCModel::SetPriorConstantAll : No parameters defined.");

   // loop over all 1-d priors
   for (int i = 0; i < nPar; ++i) {
      if (fPriorContainer[i]) {
         delete fPriorContainer[i];
         fPriorContainer[i]=0;
      }
      fPriorContainerConstant[i] = true;
   }

   RecalculatePriorConstant();

   // reset all results
   ResetResults();

   // no error
   return 1;
}

int BCModel::SetPriorGauss	(	int	index,
		double	mean,
		double	sigma
	)

Set Gaussian prior for a parameter.

Parameters:

	index	The parameter index
	mean	The mean of the Gaussian
	sigma	The sigma of the Gaussian

Returns:: An error code.

Definition at line 1738 of file BCModel.cxx.

{
   // check index range
   if (index < 0 || index >= int(GetNParameters())) {
      BCLog::OutError("BCModel::SetPriorGauss : Index out of range."); 
      return 0;
   }

   // create new function
   TF1 * f = new TF1(Form("prior_%s", GetParameter(index)->GetName().c_str()),
                     "1./sqrt(2.*TMath::Pi())/[1] * exp(- (x-[0])*(x-[0])/2./[1]/[1])",
                     GetParameter(index)->GetLowerLimit(), 
                     GetParameter(index)->GetUpperLimit()); 
   f->SetParameter(0, mean); 
   f->SetParameter(1, sigma);

   // set prior
   return SetPrior(index, f); 
}

int BCModel::SetPriorGauss	(	const char *	name,
		double	mean,
		double	sigma
	)

Set Gaussian prior for a parameter.

Parameters:

	name	The parameter name
	mean	The mean of the Gaussian
	sigma	The sigma of the Gaussian

Returns:: An error code.

Definition at line 1759 of file BCModel.cxx.

{
   // find index
   int index = -1;
   for (unsigned int i = 0; i < GetNParameters(); i++)
      if (name == GetParameter(i)->GetName())
         index = i;

   // set prior
   return SetPriorGauss(index, mean, sigma); 
}

int BCModel::SetPriorGauss	(	const char *	name,
		double	mean,
		double	sigmadown,
		double	sigmaup
	)

Set Gaussian prior for a parameter with two different widths.

Parameters:

	name	The parameter name
	mean	The mean of the Gaussian
	sigmadown	The sigma (down) of the Gaussian
	sigmaup	The sigma (up)of the Gaussian

Returns:: An error code.

Definition at line 1797 of file BCModel.cxx.

{
   // find index
   int index = -1;
   for (unsigned int i = 0; i < GetNParameters(); i++)
      if (name == GetParameter(i)->GetName())
        index = i;

   // set prior
   return SetPriorGauss(index, mean, sigmadown, sigmaup); 
}

int BCModel::SetPriorGauss	(	int	index,
		double	mean,
		double	sigmadown,
		double	sigmaup
	)

Set Gaussian prior for a parameter with two different widths.

Parameters:

	index	The parameter index
	mean	The mean of the Gaussian
	sigmadown	The sigma (down) of the Gaussian
	sigmaup	The sigma (up)of the Gaussian

Returns:: An error code.

Definition at line 1772 of file BCModel.cxx.

{
   // check index range
   if (index < 0 || index >= int(GetNParameters())) {
      BCLog::OutError("BCModel::SetPriorGauss : Index out of range."); 
      return 0;
   }

   // create new function
   TF1 * f = new TF1(Form("prior_%s", GetParameter(index)->GetName().c_str()),
                     BCMath::SplitGaussian,
                     GetParameter(index)->GetLowerLimit(),
                     GetParameter(index)->GetUpperLimit(),
                     3);
   f->SetParameter(0, mean); 
   f->SetParameter(1, sigmadown);
   f->SetParameter(2, sigmaup);

   // set prior
   return SetPrior(index, f); 

   return 0;
}

void BCModel::SetSingleDataPoint	(	BCDataSet *	dataset,
		unsigned int	index
	)

Definition at line 413 of file BCModel.cxx.

{
   if (index < 0 || index > dataset->GetNDataPoints())
      return;

   SetSingleDataPoint(dataset->GetDataPoint(index));
}

void BCModel::SetSingleDataPoint ( BCDataPoint * datapoint )

Sets a single data point as data set.

Parameters:

datapoint

A data point

Definition at line 400 of file BCModel.cxx.

{
   // create new data set consisting of a single data point
   BCDataSet * dataset = new BCDataSet();

   // add the data point
   dataset->AddDataPoint(datapoint);

   // set this new data set
   SetDataSet(dataset);
}

BCDataPoint * BCModel::VectorToDataPoint ( std::vector< double > data ) [private]

Converts a vector of doubles into a BCDataPoint

Definition at line 2277 of file BCModel.cxx.

{
   int sizeofvector = int(data.size());
   BCDataPoint * datapoint = new BCDataPoint(sizeofvector);
   datapoint->SetValues(data);
   return datapoint;
}

void BCModel::WriteMode ( const char * file )

Write mode into file

Definition at line 741 of file BCModel.cxx.

{
   ofstream ofi(file);
   if (!ofi.is_open()) {
      std::cerr << "Couldn't open file " << file << std::endl;
      return;
   }

   int npar = fParameterSet->size();
   for (int i = 0; i < npar; i++)
      ofi << fBestFitParameters.at(i) << std::endl;

   ofi << std::endl;
   ofi << "#######################################################################"
       << std::endl;
   ofi << "#" << std::endl;
   ofi << "#  This file was created automatically by BCModel::WriteMode() call."
       << std::endl;
   ofi << "#  It can be read in by call to BCModel::ReadMode()." << std::endl;
   ofi << "#  Do not modify it unless you know what you're doing." << std::endl;
   ofi << "#" << std::endl;
   ofi << "#######################################################################"
       << std::endl;
   ofi << "#" << std::endl;
   ofi << "#  Best fit parameters (mode) for model:" << std::endl;
   ofi << "#  \'" << fName.data() << "\'" << std::endl;
   ofi << "#" << std::endl;
   ofi << "#  Number of parameters: " << npar << std::endl;
   ofi << "#  Parameters ordered as above:" << std::endl;

   for (int i = 0; i < npar; i++) {
      ofi << "#     " << i << ": ";
      ofi << fParameterSet->at(i)->GetName().data() << " = ";
      ofi << fBestFitParameters.at(i) << std::endl;
   }

   ofi << "#" << std::endl;
   ofi << "########################################################################"
       << std::endl;
}