A class for handling numerical operations for models. More...

#include <BCIntegrate.h>

Inheritance diagram for BCIntegrate:

Collaboration diagram for BCIntegrate:

Public Types
Enumerators

enum	BCIntegrationMethod { kIntMonteCarlo, kIntImportance, kIntMetropolis, kIntCuba }
enum	BCMarginalizationMethod { kMargMonteCarlo, kMargMetropolis }
enum	BCOptimizationMethod { kOptSA, kOptMetropolis, kOptMinuit }
enum	BCSASchedule { kSACauchy, kSABoltzmann, kSACustom }
Public Member Functions
Constructors and destructors

	BCIntegrate ()
	BCIntegrate (BCParameterSet *par)
virtual	~BCIntegrate ()
Member functions (get)

BCIntegrate::BCIntegrationMethod	GetIntegrationMethod ()
BCIntegrate::BCMarginalizationMethod	GetMarginalizationMethod ()
BCIntegrate::BCOptimizationMethod	GetOptimizationMethod ()
BCIntegrate::BCOptimizationMethod	GetOptimizationMethodMode ()
BCIntegrate::BCSASchedule	GetSASchedule ()
void	GetRandomVector (std::vector< double > &x)
virtual void	GetRandomVectorMetro (std::vector< double > &x)
double	GetRandomPoint (std::vector< double > &x)
double	GetRandomPointImportance (std::vector< double > &x)
void	GetRandomPointMetro (std::vector< double > &x)
void	GetRandomPointSamplingMetro (std::vector< double > &x)
int	GetNiterationsPerDimension ()
int	GetNSamplesPer2DBin ()
int	GetNvar ()
int	GetNIterationsMax ()
int	GetNIterations ()
double	GetRelativePrecision ()
double	GetError ()
int	GetNbins ()
TMinuit *	GetMinuit ()
int	GetMinuitErrorFlag ()
TTree *	GetMarkovChainTree ()
std::vector< double > *	GetMarkovChainPoint ()
int *	GetMCMCIteration ()
double *	GetMarkovChainValue ()
double	GetSAT0 ()
double	GetSATmin ()
Member functions (set)

void	SetMinuitArlist (double *arglist)
void	SetFlagIgnorePrevOptimization (bool flag)
void	SetParameters (BCParameterSet *par)
void	SetVarList (int *varlist)
void	SetVar (int index)
void	SetIntegrationMethod (BCIntegrate::BCIntegrationMethod method)
void	SetMarginalizationMethod (BCIntegrate::BCMarginalizationMethod method)
void	SetOptimizationMethod (BCIntegrate::BCOptimizationMethod method)
void	SetSASchedule (BCIntegrate::BCSASchedule schedule)
void	SetNiterationsPerDimension (int niterations)
void	SetNSamplesPer2DBin (int n)
void	SetNIterationsMax (int niterations)
void	SetRelativePrecision (double relprecision)
void	SetNbins (int nbins, int index=-1)
void	SetFillErrorBand (bool flag=true)
void	UnsetFillErrorBand ()
void	SetFitFunctionIndexX (int index)
void	SetFitFunctionIndexY (int index)
void	SetFitFunctionIndices (int indexx, int indexy)
void	SetDataPointLowerBoundaries (BCDataPoint *datasetlowerboundaries)
void	SetDataPointUpperBoundaries (BCDataPoint *datasetupperboundaries)
void	SetDataPointLowerBoundary (int index, double lowerboundary)
void	SetDataPointUpperBoundary (int index, double upperboundary)
void	WriteMarkovChain (bool flag)
void	SetMarkovChainTree (TTree *tree)
void	SetMarkovChainInitialPosition (std::vector< double > position)
void	SetMarkovChainStepSize (double stepsize)
void	SetMarkovChainNIterations (int niterations)
void	SetMarkovChainAutoN (bool flag)
void	SetMode (std::vector< double > mode)
void	SetErrorBandHisto (TH2D *h)
void	SetSAT0 (double T0)
void	SetSATmin (double Tmin)
void	SetFlagWriteSAToFile (bool flag)
void	SetSATree (TTree *tree)
TTree *	GetSATree ()
void	InitializeSATree ()
Protected Attributes
std::vector< double >	fBestFitParameterErrors
std::vector< double >	fBestFitParameters
std::vector< double >	fBestFitParametersMarginalized
std::vector< bool >	fDataFixedValues
BCDataPoint *	fDataPointLowerBoundaries
BCDataPoint *	fDataPointUpperBoundaries
bool	fErrorBandContinuous
int	fErrorBandNbinsX
int	fErrorBandNbinsY
std::vector< double >	fErrorBandX
TH2D *	fErrorBandXY
bool	fFillErrorBand
int	fFitFunctionIndexX
int	fFitFunctionIndexY
double	fFlagIgnorePrevOptimization
bool	fFlagWriteMarkovChain
bool	fFlagWriteSAToFile
std::vector< TH1D * >	fHProb1D
std::vector< TH2D * >	fHProb2D
bool	fMarkovChainAutoN
int	fMarkovChainNIterations
double	fMarkovChainStepSize
TTree *	fMarkovChainTree
int	fMCMCIteration
TMinuit *	fMinuit
double	fMinuitArglist [2]
int	fMinuitErrorFlag
int	fNbins
int	fNSamplesPer2DBin
int	fNvar
double	fSALogProb
int	fSANIterations
double	fSAT0
double	fSATemperature
double	fSATmin
std::vector< double >	fSAx
TTree *	fTreeSA
Private Member Functions
void	MCMCFillErrorBand ()
void	MCMCIterationInterface ()
Private Attributes
double	fError
BCIntegrate::BCIntegrationMethod	fIntegrationMethod
BCIntegrate::BCMarginalizationMethod	fMarginalizationMethod
double	fMarkovChainValue
double *	fMax
double *	fMin
int	fNacceptedMCMC
int	fNIterations
int	fNIterationsMax
int	fNiterPerDimension
int	fNmetro
BCIntegrate::BCOptimizationMethod	fOptimizationMethod
BCIntegrate::BCOptimizationMethod	fOptimizationMethodMode
double	fRelativePrecision
BCIntegrate::BCSASchedule	fSASchedule
int *	fVarlist
BCParameterSet *	fx
std::vector< double >	fXmetro0
std::vector< double >	fXmetro1
Member functions (miscellaneous methods)

void	DeleteVarList ()
void	ResetVarlist (int v)
void	UnsetVar (int index)
virtual double	Eval (std::vector< double > x)
virtual double	LogEval (std::vector< double > x)
virtual double	EvalSampling (std::vector< double > x)
double	LogEvalSampling (std::vector< double > x)
double	EvalPrint (std::vector< double > x)
virtual double	FitFunction (std::vector< double >, std::vector< double >)
double	Integrate ()
double	IntegralMC (std::vector< double > x, int *varlist)
double	IntegralMC (std::vector< double > x)
double	IntegralMetro (std::vector< double > x)
double	IntegralImportance (std::vector< double > x)
double	CubaIntegrate (int method, std::vector< double > parameters_double, std::vector< double > parameters_int)
double	CubaIntegrate ()
TH1D *	Marginalize (BCParameter *parameter)
TH2D *	Marginalize (BCParameter parameter1, BCParameter parameter2)
TH1D *	MarginalizeByIntegrate (BCParameter *parameter)
TH2D *	MarginalizeByIntegrate (BCParameter parameter1, BCParameter parameter2)
TH1D *	MarginalizeByMetro (BCParameter *parameter)
TH2D *	MarginalizeByMetro (BCParameter parameter1, BCParameter parameter2)
int	MarginalizeAllByMetro (const char *name)
TH1D *	GetH1D (int parIndex)
int	GetH2DIndex (int parIndex1, int parIndex2)
TH2D *	GetH2D (int parIndex1, int parIndex2)
void	InitMetro ()
void	SAInitialize ()
virtual void	FindModeMinuit (std::vector< double > start=std::vector< double >(0), int printlevel=1)
void	FindModeMCMC ()
void	FindModeSA (std::vector< double > start=std::vector< double >(0))
double	SATemperature (double t)
double	SATemperatureBoltzmann (double t)
double	SATemperatureCauchy (double t)
virtual double	SATemperatureCustom (double t)
std::vector< double >	GetProposalPointSA (std::vector< double > x, int t)
std::vector< double >	GetProposalPointSABoltzmann (std::vector< double > x, int t)
std::vector< double >	GetProposalPointSACauchy (std::vector< double > x, int t)
virtual std::vector< double >	GetProposalPointSACustom (std::vector< double > x, int t)
std::vector< double >	SAHelperGetRandomPointOnHypersphere ()
double	SAHelperGetRadialCauchy ()
double	SAHelperSinusToNIntegral (int dim, double theta)
virtual void	MCMCUserIterationInterface ()
int	IntegrateResetResults ()
static void	CubaIntegrand (const int ndim, const double xx[], const int ncomp, double ff[])
static void	FCNLikelihood (int &npar, double grad, double &fval, double par, int flag)

Detailed Description

A class for handling numerical operations for models.

Author:: Daniel Kollar; Kevin Kröninger

Version:: 1.0

Date:: 08.2008 This is a base class for a model class. It contains numerical methods to carry out the integration, marginalization, peak finding etc.

Definition at line 45 of file BCIntegrate.h.

Member Enumeration Documentation

enum BCIntegrate::BCIntegrationMethod

An enumerator for the integration algorithm

Enumerator:

kIntMonteCarlo
kIntImportance
kIntMetropolis
kIntCuba

Definition at line 55 of file BCIntegrate.h.

{ kIntMonteCarlo, kIntImportance, kIntMetropolis, kIntCuba };

enum BCIntegrate::BCMarginalizationMethod

An enumerator for the marginalization algorithm

Enumerator:

kMargMonteCarlo
kMargMetropolis

Definition at line 59 of file BCIntegrate.h.

{ kMargMonteCarlo, kMargMetropolis };

enum BCIntegrate::BCOptimizationMethod

An enumerator for the mode finding algorithm

Enumerator:

kOptSA
kOptMetropolis
kOptMinuit

Definition at line 63 of file BCIntegrate.h.

{ kOptSA, kOptMetropolis, kOptMinuit };

enum BCIntegrate::BCSASchedule

An enumerator for the Simulated Annealing schedule

Enumerator:

kSACauchy
kSABoltzmann
kSACustom

Definition at line 67 of file BCIntegrate.h.

{ kSACauchy, kSABoltzmann, kSACustom };

Constructor & Destructor Documentation

BCIntegrate::BCIntegrate ( )

The default constructor

Definition at line 31 of file BCIntegrate.cxx.

                         : BCEngineMCMC()
{
   fNvar=0;
   fNiterPerDimension = 100;
   fNSamplesPer2DBin = 100;

   fNIterationsMax    = 1000000;
   fRelativePrecision = 1e-3;

#ifdef HAVE_CUBA_H
   fIntegrationMethod = BCIntegrate::kIntCuba;
#else
   fIntegrationMethod = BCIntegrate::kIntMonteCarlo;
#endif
   fMarginalizationMethod = BCIntegrate::kMargMetropolis;
   fOptimizationMethod = BCIntegrate::kOptMinuit;
   fOptimizationMethodMode = BCIntegrate::kOptMinuit;

   fNbins=100;

   fDataPointLowerBoundaries = 0;
   fDataPointUpperBoundaries = 0;

   fFillErrorBand = false;

   fFitFunctionIndexX = -1;
   fFitFunctionIndexY = -1;

   fErrorBandNbinsX = 100;
   fErrorBandNbinsY = 500;

   fErrorBandContinuous = true;

   fMinuit = 0;
   fMinuitArglist[0] = 20000;
   fMinuitArglist[1] = 0.01;
   fFlagIgnorePrevOptimization = false;

   fFlagWriteMarkovChain = false;
   fMarkovChainTree = 0;

   fMarkovChainStepSize = 0.1;

   fMarkovChainAutoN = true;

   fSAT0 = 100.0;
   fSATmin = 0.1;
   fSASchedule = BCIntegrate::kSACauchy;
   fFlagWriteSAToFile = false;
}

BCIntegrate::BCIntegrate ( BCParameterSet * par )

A constructor

Definition at line 83 of file BCIntegrate.cxx.

                                             : BCEngineMCMC(1)
{
   fNvar=0;
   fNiterPerDimension = 100;
   fNSamplesPer2DBin = 100;

   fNbins=100;

   SetParameters(par);

   fDataPointLowerBoundaries = 0;
   fDataPointUpperBoundaries = 0;

   fFillErrorBand = false;

   fFitFunctionIndexX = -1;
   fFitFunctionIndexY = -1;

   fErrorBandNbinsX = 100;
   fErrorBandNbinsY = 200;

   fErrorBandContinuous = true;

   fMinuit = 0;
   fMinuitArglist[0] = 20000;
   fMinuitArglist[1] = 0.01;
   fFlagIgnorePrevOptimization = false;

   fFlagWriteMarkovChain = false;
   fMarkovChainTree = 0;

   fMarkovChainStepSize = 0.1;

   fMarkovChainAutoN = true;

   fSAT0 = 100.0;
   fSATmin = 0.1;
   fTreeSA = false;
}

BCIntegrate::~BCIntegrate ( ) [virtual]

The default destructor

Definition at line 124 of file BCIntegrate.cxx.

{
   DeleteVarList();

   fx=0;

   if (fMinuit)
      delete fMinuit;

   int n1 = fHProb1D.size();
   if(n1>0) {
      for (int i=0;i<n1;i++)
         delete fHProb1D.at(i);
   }

   int n2 = fHProb2D.size();
   if(n2>0) {
      for (int i=0;i<n2;i++)
         delete fHProb2D.at(i);
   }
}

Member Function Documentation

void BCIntegrate::CubaIntegrand	(	const int *	ndim,
		const double	xx[],
		const int *	ncomp,
		double	ff[]
	)			`[static]`

Integrand for the Cuba library.

Parameters:

	ndim	The number of dimensions to integrate over
	xx	The point in parameter space to integrate over (scaled to 0 - 1 per dimension)
	ncomp	The number of components of the integrand (usually 1)
	ff	The function value

Returns:: The integral

Definition at line 1775 of file BCIntegrate.cxx.

{
#ifdef HAVE_CUBA_H
   // scale variables
   double jacobian = 1.0;

   std::vector<double> scaled_parameters;

   for (int i = 0; i < *ndim; i++) {
      double range = global_this->fx->at(i)->GetUpperLimit() -  global_this->fx->at(i)->GetLowerLimit();

      // multiply range to jacobian
      jacobian *= range;

      // get the scaled parameter value
      scaled_parameters.push_back(global_this->fx->at(i)->GetLowerLimit() + xx[i] * range);
   }

   // call function to integrate
   ff[0] = global_this->Eval(scaled_parameters);

   // multiply jacobian
   ff[0] *= jacobian;

   // multiply fudge factor
   ff[0] *= 1e99;

   // remove parameter vector
   scaled_parameters.clear();
#else
   BCLog::OutError("!!! This version of BAT is compiled without Cuba.");
   BCLog::OutError("    Use other integration methods or install Cuba and recompile BAT.");
#endif
}

double BCIntegrate::CubaIntegrate	(	int	method,
		std::vector< double >	parameters_double,
		std::vector< double >	parameters_int
	)

Calculate integral using the Cuba library. For details see documentation.

Parameters:

	method	A short cut for the method
	parameters_double	A vector of parameters (double)
	parameters_int	A vector of parameters (int)

Returns:: The integral

Definition at line 1849 of file BCIntegrate.cxx.

{
#ifdef HAVE_CUBA_H
   const int NDIM      = int(fx ->size());
   const int NCOMP     = 1;

   const double EPSREL = parameters_double[0];
   const double EPSABS = parameters_double[1];
   const int VERBOSE   = int(parameters_int[0]);
   const int MINEVAL   = int(parameters_int[1]);
   const int MAXEVAL   = int(parameters_int[2]);

   int neval;
   int fail;
   int nregions;
   double integral[NCOMP];
   double error[NCOMP];
   double prob[NCOMP];

   global_this = this;

   integrand_t an_integrand = &BCIntegrate::CubaIntegrand;

   if (method == 0) {
      // set VEGAS specific parameters
      const int NSTART    = int(parameters_int[3]);
      const int NINCREASE = int(parameters_int[4]);

      // call VEGAS integration method
      Vegas(NDIM, NCOMP, an_integrand,
         EPSREL, EPSABS, VERBOSE, MINEVAL, MAXEVAL,
         NSTART, NINCREASE,
         &neval, &fail, integral, error, prob);
   }
   else if (method == 1) {
      // set SUAVE specific parameters
      const int LAST     = int(parameters_int[3]);
      const int NNEW     = int(parameters_int[4]);
      const int FLATNESS = int(parameters_int[5]);

      // call SUAVE integration method
      Suave(NDIM, NCOMP, an_integrand,
         EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
         NNEW, FLATNESS,
         &nregions, &neval, &fail, integral, error, prob);
   }
   else if (method == 2) {
      // set DIVONNE specific parameters
      const int KEY1         = int(parameters_int[3]);
      const int KEY2         = int(parameters_int[4]);
      const int KEY3         = int(parameters_int[5]);
      const int MAXPASS      = int(parameters_int[6]);
      const int BORDER       = int(parameters_int[7]);
      const int MAXCHISQ     = int(parameters_int[8]);
      const int MINDEVIATION = int(parameters_int[9]);
      const int NGIVEN       = int(parameters_int[10]);
      const int LDXGIVEN     = int(parameters_int[11]);
      const int NEXTRA       = int(parameters_int[12]);

      // call DIVONNE integration method
      Divonne(NDIM, NCOMP, an_integrand,
         EPSREL, EPSABS, VERBOSE, MINEVAL, MAXEVAL,
         KEY1, KEY2, KEY3, MAXPASS, BORDER, MAXCHISQ, MINDEVIATION,
         NGIVEN, LDXGIVEN, NULL, NEXTRA, NULL,
         &nregions, &neval, &fail, integral, error, prob);
   }
   else if (method == 3) {
      // set CUHRE specific parameters
      const int LAST = int(parameters_int[3]);
      const int KEY  = int(parameters_int[4]);

      // call CUHRE integration method
      Cuhre(NDIM, NCOMP, an_integrand,
         EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL, KEY,
         &nregions, &neval, &fail, integral, error, prob);
   }
   else {
      std::cout << " Integration method not available. " << std::endl;
      integral[0] = -1e99;
   }

   if (fail != 0) {
      std::cout << " Warning, integral did not converge with the given set of parameters. "<< std::endl;
      std::cout << " neval    = " << neval       << std::endl;
      std::cout << " fail     = " << fail        << std::endl;
      std::cout << " integral = " << integral[0] << std::endl;
      std::cout << " error    = " << error[0]    << std::endl;
      std::cout << " prob     = " << prob[0]     << std::endl;
   }

   return integral[0] / 1e99;
#else
   BCLog::OutError("!!! This version of BAT is compiled without Cuba.");
   BCLog::OutError("    Use other integration methods or install Cuba and recompile BAT.");
   return 0.;
#endif
}

double BCIntegrate::CubaIntegrate ( )

Definition at line 1812 of file BCIntegrate.cxx.

{
#ifdef HAVE_CUBA_H
   double EPSREL = 1e-3;
   double EPSABS = 1e-12;
   double VERBOSE   = 0;
   double MINEVAL   = 0;
   double MAXEVAL   = 2000000;
   double NSTART    = 25000;
   double NINCREASE = 25000;

   std::vector<double> parameters_double;
   std::vector<double>    parameters_int;

   parameters_double.push_back(EPSREL);
   parameters_double.push_back(EPSABS);

   parameters_int.push_back(VERBOSE);
   parameters_int.push_back(MINEVAL);
   parameters_int.push_back(MAXEVAL);
   parameters_int.push_back(NSTART);
   parameters_int.push_back(NINCREASE);

   // print to log
   BCLog::OutDebug( Form(" --> Running Cuba/Vegas integation over %i dimensions.", fNvar));
   BCLog::OutDebug( Form(" --> Maximum number of iterations: %i", (int)MAXEVAL));
   BCLog::OutDebug( Form(" --> Aimed relative precision:     %e", EPSREL));

   return CubaIntegrate(0, parameters_double, parameters_int);
#else
   BCLog::OutError("!!! This version of BAT is compiled without Cuba.");
   BCLog::OutError("    Use other integration methods or install Cuba and recompile BAT.");
   return 0.;
#endif
}

void BCIntegrate::DeleteVarList ( )

Frees the memory for integration variables

Definition at line 198 of file BCIntegrate.cxx.

{
   if(fNvar) {
      delete[] fVarlist;
      fVarlist=0;

      delete[] fMin;
      fMin=0;

      delete[] fMax;
      fMax=0;

      fx=0;
      fNvar=0;
   }
}

double BCIntegrate::Eval ( std::vector< double > x ) [virtual]

Evaluate the unnormalized probability at a point in parameter space. Method needs to be overloaded by the user.

Parameters:

x

The point in parameter space

Returns:: The unnormalized probability

Reimplemented in BCModel.

Definition at line 291 of file BCIntegrate.cxx.

{
   BCLog::OutWarning( "BCIntegrate::Eval. No function. Function needs to be overloaded.");
   return 0;
}

double BCIntegrate::EvalPrint ( std::vector< double > x )

Evaluate the un-normalized probability at a point in parameter space and prints the result to the log.

Parameters:

x

The point in parameter space

Returns:: The un-normalized probability

See also:: Eval(std::vector <double> x)

Definition at line 318 of file BCIntegrate.cxx.

{
   double val=Eval(x);
   BCLog::OutDebug(Form("BCIntegrate::EvalPrint. Value: %d.", val));

   return val;
}

double BCIntegrate::EvalSampling ( std::vector< double > x ) [virtual]

Evaluate the sampling function at a point in parameter space. Method needs to be overloaded by the user.

Parameters:

x

The point in parameter space

Returns:: The value of the sampling function

Reimplemented in BCModel.

Definition at line 305 of file BCIntegrate.cxx.

{
   BCLog::OutWarning( "BCIntegrate::EvalSampling. No function. Function needs to be overloaded.");
   return 0;
}

void BCIntegrate::FCNLikelihood	(	int &	npar,
		double *	grad,
		double &	fval,
		double *	par,
		int	flag
	)			`[static]`

Definition at line 1718 of file BCIntegrate.cxx.

{
   // copy parameters
   std::vector <double> parameters;

   int n = global_this->GetNvar();

   for (int i = 0; i < n; i++)
      parameters.push_back(par[i]);

   fval = - global_this->LogEval(parameters);

   // delete parameters
   parameters.clear();
}

void BCIntegrate::FindModeMCMC ( )

Does the mode finding using Markov Chain Monte Carlo (prerun only!)

Definition at line 1737 of file BCIntegrate.cxx.

{
   // call PreRun
//   if (flag_run == 0)
   if (!fMCMCFlagPreRun)
      MCMCMetropolisPreRun();

   // find global maximum
   //   double probmax = (MCMCGetMaximumLogProb()).at(0);

   double probmax = 0;

   if ( int(fBestFitParameters.size()) == fNvar) {
      probmax = Eval( fBestFitParameters );
//      fBestFitParameters = MCMCGetMaximumPoint(0);
   }

   // loop over all chains and find the maximum point
   for (int i = 0; i < fMCMCNChains; ++i) {
      double prob = exp( (MCMCGetMaximumLogProb()).at(i));

      // copy the point into the vector
      if ( (prob >= probmax && !fFlagIgnorePrevOptimization) || fFlagIgnorePrevOptimization) {
         probmax = prob;

         fBestFitParameters.clear();
         fBestFitParameterErrors.clear();
         fBestFitParameters = MCMCGetMaximumPoint(i);

         for (int j = 0; j < fNvar; ++j)
            fBestFitParameterErrors.push_back(-1.); // error undefined

         fOptimizationMethodMode = BCIntegrate::kOptMetropolis;
      }
   }
}

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

Does the mode finding 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 in BCModel.

Definition at line 1234 of file BCIntegrate.cxx.

{
   bool have_start = true;

   // check start values
   if (int(start.size()) != fNvar)
      have_start = false;

   // set global this
   global_this = this;

   // define minuit
   if (fMinuit)
      delete fMinuit;
   fMinuit = new TMinuit(fNvar);

   // set print level
   fMinuit->SetPrintLevel(printlevel);

   // set function
   fMinuit->SetFCN(&BCIntegrate::FCNLikelihood);

   // set UP for likelihood
   fMinuit->SetErrorDef(0.5);

   // set parameters
   int flag;
   for (int i = 0; i < fNvar; i++) {
      double starting_point = (fMin[i]+fMax[i])/2.;
      if(have_start)
         starting_point = start[i];
         fMinuit->mnparm(i,
                         fx->at(i)->GetName().data(),
                         starting_point,
                         (fMax[i]-fMin[i])/100.,
                         fMin[i],
                         fMax[i],
                         flag);
   }

   // do mcmc minimization
   //   fMinuit->mnseek();

   // do minimization
   fMinuit->mnexcm("MIGRAD", fMinuitArglist, 2, flag);

   // improve search for local minimum
   //   fMinuit->mnimpr();

   // copy flag
   fMinuitErrorFlag = flag;

   // check if mode found by minuit is better than previous estimation
   double probmax = 0;
   bool valid = false;

   if ( int(fBestFitParameters.size()) == fNvar) {
      valid = true;
      for (int i = 0; i < fNvar; ++i)
         if (fBestFitParameters.at(i) < fMin[i] || fBestFitParameters.at(i) > fMax[i])
            valid= false;

      if (valid)
         probmax = Eval( fBestFitParameters );
   }

   std::vector<double> tempvec;
   for (int i = 0; i < fNvar; i++) {
      double par;
      double parerr;
      fMinuit->GetParameter(i, par, parerr);
      tempvec.push_back(par);
   }
   double probmaxminuit = Eval( tempvec );

   // set best fit parameters
   if ( (probmaxminuit > probmax && !fFlagIgnorePrevOptimization) || !valid || fFlagIgnorePrevOptimization) {
      fBestFitParameters.clear();
      fBestFitParameterErrors.clear();

      for (int i = 0; i < fNvar; i++) {
         double par;
         double parerr;
         fMinuit->GetParameter(i, par, parerr);
         fBestFitParameters.push_back(par);
         fBestFitParameterErrors.push_back(parerr);
      }

      // set optimization method used to find the mode
      fOptimizationMethodMode = BCIntegrate::kOptMinuit;
   }

   // delete minuit
   delete fMinuit;
   fMinuit = 0;

   return;
}

void BCIntegrate::FindModeSA ( std::vector< double > start = std::vector<double>(0) )

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

Parameters:

start

point in parameter space from thich the mode finding is started.

Definition at line 1349 of file BCIntegrate.cxx.

{
   // note: if f(x) is the function to be minimized, then
   // f(x) := - LogEval(parameters)

   bool have_start = true;
   std::vector<double> x, y, best_fit; // vectors for current state, new proposed state and best fit up to now
   double fval_x, fval_y, fval_best_fit; // function values at points x, y and best_fit (we save them rather than to re-calculate them every time)
   int t = 1; // time iterator

   // check start values
   if (int(start.size()) != fNvar)
      have_start = false;

   // if no starting point is given, set to center of parameter space
   if ( !have_start ) {
      start.clear();
      for (int i = 0; i < fNvar; i++)
         start.push_back((fMin[i]+fMax[i])/2.);
   }

   // set current state and best fit to starting point
   x.clear();
   best_fit.clear();
   for (int i = 0; i < fNvar; i++) {
      x.push_back(start[i]);
      best_fit.push_back(start[i]);
   }
   // calculate function value at starting point
   fval_x = fval_best_fit = LogEval(x);

   // run while still "hot enough"
   while ( SATemperature(t) > fSATmin ) {
      // generate new state
      y = GetProposalPointSA(x, t);

      // check if the proposed point is inside the phase space
      // if not, reject it
      bool is_in_ranges = true;

      for (int i = 0; i < fNvar; i++)
         if (y[i] > fMax[i] || y[i] < fMin[i])
            is_in_ranges = false;

      if ( !is_in_ranges )
         ; // do nothing...
      else {
         // calculate function value at new point
         fval_y = LogEval(y);

         // is it better than the last one?
         // if so, update state and chef if it is the new best fit...
         if (fval_y >= fval_x) {
            x.clear();
            for (int i = 0; i < fNvar; i++)
               x.push_back(y[i]);

            fval_x = fval_y;

            if (fval_y > fval_best_fit) {
               best_fit.clear();
               for (int i = 0; i < fNvar; i++)
                  best_fit.push_back(y[i]);

               fval_best_fit = fval_y;
            }
         }
         // ...else, only accept new state w/ certain probability
         else {
            if (fRandom->Rndm() <= exp( (fval_y - fval_x) / SATemperature(t) )) {
               x.clear();
               for (int i = 0; i < fNvar; i++)
                  x.push_back(y[i]);

               fval_x = fval_y;
            }
         }
      }

      // update tree variables
      fSANIterations = t;
      fSATemperature = SATemperature(t);
      fSALogProb = fval_x;
      fSAx.clear();
      for (int i = 0; i < fNvar; ++i)
         fSAx.push_back(x[i]);

      // fill tree
      if (fFlagWriteSAToFile)
         fTreeSA->Fill();

      // increate t
      t++;
   }

   // check if mode found by minuit is better than previous estimation
   double probmax = 0;
   bool valid = false;

   if ( int(fBestFitParameters.size()) == fNvar) {
      valid = true;
      for (int i = 0; i < fNvar; ++i)
         if (fBestFitParameters.at(i) < fMin[i] || fBestFitParameters.at(i) > fMax[i])
            valid= false;

      if (valid)
         probmax = Eval( fBestFitParameters );
   }

   double probmaxsa = Eval( best_fit);

   if ( (probmaxsa > probmax  && !fFlagIgnorePrevOptimization) || !valid || fFlagIgnorePrevOptimization) {
      // set best fit parameters
      fBestFitParameters.clear();
      fBestFitParameterErrors.clear();

      for (int i = 0; i < fNvar; i++) {
         fBestFitParameters.push_back(best_fit[i]);
         fBestFitParameterErrors.push_back(-1.); // error undefined
      }

      // set optimization moethod used to find the mode
      fOptimizationMethodMode = BCIntegrate::kOptSA;
   }

   return;
}

virtual double BCIntegrate::FitFunction	(	std::vector< double >	,
		std::vector< double >
	)			`[inline, virtual]`

Defines a fit function.

Parameters:

	parameters	A set of parameter values
	x	A vector of x-values

Returns:: The value of the fit function at the x-values given a set of parameters

Reimplemented in BCEfficiencyFitter, BCGraphFitter, and BCHistogramFitter.

Definition at line 489 of file BCIntegrate.h.

         { return 0.; };

double BCIntegrate::GetError ( ) [inline]

Definition at line 181 of file BCIntegrate.h.

         { return fError; };

TH1D * BCIntegrate::GetH1D ( int parIndex )

Parameters:

parIndex1

Index of parameter

Returns:: Pointer to 1D histogram (TH1D) of marginalized distribution wrt. parameter with given index.

Definition at line 968 of file BCIntegrate.cxx.

{
   if(fHProb1D.size()==0) {
      BCLog::OutWarning("BCModel::GetH1D. MarginalizeAll() has to be run prior to this to fill the distributions.");
      return 0;
   }

   if(parIndex<0 || parIndex>fNvar) {
      BCLog::OutWarning(Form("BCIntegrate::GetH1D. Parameter index %d is invalid.",parIndex));
      return 0;
   }

   return fHProb1D[parIndex];
}

TH2D * BCIntegrate::GetH2D	(	int	parIndex1,
		int	parIndex2
	)

Parameters:

	parIndex1	Index of first parameter
	parIndex2	Index of second parameter, with parIndex2>parIndex1

Returns:: Pointer to 2D histogram (TH2D) of marginalized distribution wrt. parameters with given indeces.

Definition at line 1020 of file BCIntegrate.cxx.

{
   if(fHProb2D.size()==0) {
      BCLog::OutWarning("BCModel::GetH2D. MarginalizeAll() has to be run prior to this to fill the distributions.");
      return 0;
   }

   int hindex = GetH2DIndex(parIndex1, parIndex2);
   if(hindex==-1)
      return 0;

   if(hindex>(int)fHProb2D.size()-1) {
      BCLog::OutWarning("BCIntegrate::GetH2D. Got invalid index from GetH2DIndex(). Something went wrong.");
      return 0;
   }

   return fHProb2D[hindex];
}

int BCIntegrate::GetH2DIndex	(	int	parIndex1,
		int	parIndex2
	)

Parameters:

	parIndex1	Index of first parameter
	parIndex2	Index of second parameter, with parIndex2>parIndex1

Returns:: Index of the distribution in the vector of 2D distributions, which corresponds to the combination of parameters with given indeces

Definition at line 984 of file BCIntegrate.cxx.

{
   if(parIndex1>fNvar-1 || parIndex1<0) {
      BCLog::OutWarning(Form("BCIntegrate::GetH2DIndex. Parameter index %d is invalid", parIndex1));
      return -1;
   }

   if(parIndex2>fNvar-1 || parIndex2<0) {
      BCLog::OutWarning(Form("BCIntegrate::GetH2DIndex. Parameter index %d is invalid", parIndex2));
      return -1;
   }

   if(parIndex1==parIndex2) {
      BCLog::OutWarning(Form("BCIntegrate::GetH2DIndex. Parameters have equal indeces: %d , %d", parIndex1, parIndex2));
      return -1;
   }

   if(parIndex1>parIndex2) {
      BCLog::OutWarning("BCIntegrate::GetH2DIndex. First parameters must be smaller than second (sorry).");
      return -1;
   }

   int index=0;
   for(int i=0;i<fNvar-1;i++)
      for(int j=i+1;j<fNvar;j++) {
         if(i==parIndex1 && j==parIndex2)
            return index;
         index++;
      }

   BCLog::OutWarning("BCIntegrate::GetH2DIndex. Invalid index combination.");

   return -1;
}

BCIntegrate::BCIntegrationMethod BCIntegrate::GetIntegrationMethod ( ) [inline]

Returns:: The integration method

Definition at line 93 of file BCIntegrate.h.

         { return fIntegrationMethod; };

BCIntegrate::BCMarginalizationMethod BCIntegrate::GetMarginalizationMethod ( ) [inline]

Returns:: The marginalization method

Definition at line 98 of file BCIntegrate.h.

         { return fMarginalizationMethod; };

std::vector<double>* BCIntegrate::GetMarkovChainPoint ( ) [inline]

Definition at line 203 of file BCIntegrate.h.

         { return &fXmetro1; };

TTree* BCIntegrate::GetMarkovChainTree ( ) [inline]

Definition at line 198 of file BCIntegrate.h.

         { return fMarkovChainTree; };

double* BCIntegrate::GetMarkovChainValue ( ) [inline]

Returns the value of the loglikelihood at the point fXmetro1

Definition at line 213 of file BCIntegrate.h.

         { return &fMarkovChainValue; };

int* BCIntegrate::GetMCMCIteration ( ) [inline]

Returns the iteration of the MCMC

Definition at line 208 of file BCIntegrate.h.

         { return &fMCMCIteration; };

TMinuit * BCIntegrate::GetMinuit ( )

Definition at line 1224 of file BCIntegrate.cxx.

{
   if (!fMinuit)
      fMinuit = new TMinuit();

   return fMinuit;
}

int BCIntegrate::GetMinuitErrorFlag ( ) [inline]

Definition at line 193 of file BCIntegrate.h.

         { return fMinuitErrorFlag; };

int BCIntegrate::GetNbins ( ) [inline]

Definition at line 186 of file BCIntegrate.h.

         { return fNbins; };

int BCIntegrate::GetNIterations ( ) [inline]

Returns:: The number of iterations for the most recent Monte Carlo integration

Definition at line 171 of file BCIntegrate.h.

         { return fNIterations; };

int BCIntegrate::GetNIterationsMax ( ) [inline]

Returns:: The number of maximum iterations for Monte Carlo integration

Definition at line 166 of file BCIntegrate.h.

         { return fNIterationsMax; };

int BCIntegrate::GetNiterationsPerDimension ( ) [inline]

Returns:: The number of iterations per dimension for the Monte Carlo integration

Definition at line 151 of file BCIntegrate.h.

         { return fNiterPerDimension; };

int BCIntegrate::GetNSamplesPer2DBin ( ) [inline]

Returns:: Number of samples per 2D bin per variable in the Metropolis marginalization.

Definition at line 156 of file BCIntegrate.h.

         { return fNSamplesPer2DBin; };

int BCIntegrate::GetNvar ( ) [inline]

Returns:: The number of variables to integrate over

Definition at line 161 of file BCIntegrate.h.

         { return fNvar; };

BCIntegrate::BCOptimizationMethod BCIntegrate::GetOptimizationMethod ( ) [inline]

Returns:: The current optimization method

Definition at line 103 of file BCIntegrate.h.

         { return fOptimizationMethod; };

BCIntegrate::BCOptimizationMethod BCIntegrate::GetOptimizationMethodMode ( ) [inline]

Returns:: The optimization method used to find the mode

Definition at line 108 of file BCIntegrate.h.

         { return fOptimizationMethodMode; };

std::vector< double > BCIntegrate::GetProposalPointSA	(	std::vector< double >	x,
		int	t
	)

Generates a new state in a neighbourhood around x that is to be accepted or rejected by the Simulated Annealing algorithm. Delegates the generation to the appropriate method according to fSASchedule.

Parameters:

	x	last state.
	t	time iterator to determine current temperature.

Definition at line 1514 of file BCIntegrate.cxx.

{
   // do we have Cauchy (default), Boltzmann or custom annealing schedule?
   if (fSASchedule == BCIntegrate::kSABoltzmann)
      return GetProposalPointSABoltzmann(x, t);
   else if (fSASchedule == BCIntegrate::kSACauchy)
      return GetProposalPointSACauchy(x, t);
   else
      return GetProposalPointSACustom(x, t);
}

std::vector< double > BCIntegrate::GetProposalPointSABoltzmann	(	std::vector< double >	x,
		int	t
	)

Generates a new state in a neighbourhood around x that is to be accepted or rejected by the Simulated Annealing algorithm. This method is used for Boltzmann annealing schedule.

Parameters:

	x	last state.
	t	time iterator to determine current temperature.

Definition at line 1527 of file BCIntegrate.cxx.

{
   std::vector<double> y;
   y.clear();
   double new_val, norm;

   for (int i = 0; i < fNvar; i++) {
      norm = (fMax[i] - fMin[i]) * SATemperature(t) / 2.;
      new_val = x[i] + norm * fRandom->Gaus();
      y.push_back(new_val);
   }

   return y;
}

std::vector< double > BCIntegrate::GetProposalPointSACauchy	(	std::vector< double >	x,
		int	t
	)

Generates a new state in a neighbourhood around x that is to be accepted or rejected by the Simulated Annealing algorithm. This method is used for Cauchy annealing schedule.

Parameters:

	x	last state.
	t	time iterator to determine current temperature.

Definition at line 1544 of file BCIntegrate.cxx.

{
   std::vector<double> y;
   y.clear();

   if (fNvar == 1) {
      double cauchy, new_val, norm;

      norm = (fMax[0] - fMin[0]) * SATemperature(t) / 2.;
      cauchy = tan(3.14159 * (fRandom->Rndm() - 0.5));
      new_val = x[0] + norm * cauchy;
      y.push_back(new_val);
   }
   else {
      // use sampling to get radial n-dim Cauchy distribution

      // first generate a random point uniformly distributed on a
      // fNvar-dimensional hypersphere
      y = SAHelperGetRandomPointOnHypersphere();

      // scale the vector by a random factor determined by the radial
      // part of the fNvar-dimensional Cauchy distribution
      double radial = SATemperature(t) * SAHelperGetRadialCauchy();

      // scale y by radial part and the size of dimension i in phase space
      // afterwards, move by x
      for (int i = 0; i < fNvar; i++)
         y[i] = (fMax[i] - fMin[i]) * y[i] * radial / 2. + x[i];
   }

   return y;
}

std::vector< double > BCIntegrate::GetProposalPointSACustom	(	std::vector< double >	x,
		int	t
	)			`[virtual]`

Generates a new state in a neighbourhood around x that is to be accepted or rejected by the Simulated Annealing algorithm. This is a virtual method to be overridden by a user-defined custom proposal function.

Parameters:

	x	last state.
	t	time iterator to determine current temperature.

Definition at line 1579 of file BCIntegrate.cxx.

{
   BCLog::OutError("BCIntegrate::GetProposalPointSACustom : No custom proposal function defined");
   return std::vector<double>(fNvar);
}

double BCIntegrate::GetRandomPoint ( std::vector< double > & x )

Fills a vector of (flat) random numbers in the limits of the parameters and returns the probability at that point

Parameters:

x

A vector of doubles

Returns:: The (unnormalized) probability at the random point

Definition at line 1040 of file BCIntegrate.cxx.

{
   GetRandomVector(x);

   for(int i=0;i<fNvar;i++)
      x[i]=fMin[i]+x[i]*(fMax[i]-fMin[i]);

   return Eval(x);
}

double BCIntegrate::GetRandomPointImportance ( std::vector< double > & x )

Fills a vector of random numbers in the limits of the parameters sampled by the sampling function and returns the probability at that point

Parameters:

x

A vector of doubles

Returns:: The (unnormalized) probability at the random point

Definition at line 1051 of file BCIntegrate.cxx.

{
   double p = 1.1;
   double g = 1.0;

   while (p>g) {
      GetRandomVector(x);

      for(int i=0;i<fNvar;i++)
         x[i] = fMin[i] + x[i] * (fMax[i]-fMin[i]);

      p = fRandom->Rndm();

      g = EvalSampling(x);
   }

   return Eval(x);
}

void BCIntegrate::GetRandomPointMetro ( std::vector< double > & x )

Fills a vector of random numbers in the limits of the parameters sampled by the probality function and returns the probability at that point (Metropolis)

Parameters:

x

A vector of doubles

Definition at line 1094 of file BCIntegrate.cxx.

{
   // get new point
   GetRandomVectorMetro(fXmetro1);

   // scale the point to the allowed region and stepsize
   int in=1;
   for(int i=0;i<fNvar;i++) {
      fXmetro1[i] = fXmetro0[i] + fXmetro1[i] * (fMax[i]-fMin[i]);

      // check whether the generated point is inside the allowed region
      if( fXmetro1[i]<fMin[i] || fXmetro1[i]>fMax[i] )
         in=0;
   }

   // calculate the log probabilities and compare old and new point

   double p0 = fMarkovChainValue; // old point
   double p1 = 0; // new point
   int accept=0;

   if(in) {
      p1 = LogEval(fXmetro1);

      if(p1>=p0)
         accept=1;
      else {
         double r=log(fRandom->Rndm());
         if(r<p1-p0)
            accept=1;
      }
   }

   // fill the return point after the decision
   if(accept) {
      // increase counter
      fNacceptedMCMC++;

      for(int i=0;i<fNvar;i++) {
         fXmetro0[i]=fXmetro1[i];
         x[i]=fXmetro1[i];
         fMarkovChainValue = p1;
      }
   }
   else
      for(int i=0;i<fNvar;i++) {
         x[i]=fXmetro0[i];
         fXmetro1[i] = x[i];
         fMarkovChainValue = p0;
      }

   fNmetro++;
}

void BCIntegrate::GetRandomPointSamplingMetro ( std::vector< double > & x )

Fills a vector of random numbers in the limits of the parameters sampled by the sampling function and returns the probability at that point (Metropolis)

Parameters:

x

A vector of doubles

Definition at line 1149 of file BCIntegrate.cxx.

{
   // get new point
   GetRandomVectorMetro(fXmetro1);

   // scale the point to the allowed region and stepsize
   int in=1;
   for(int i=0;i<fNvar;i++) {
      fXmetro1[i] = fXmetro0[i] + fXmetro1[i] * (fMax[i]-fMin[i]);

      // check whether the generated point is inside the allowed region
      if( fXmetro1[i]<fMin[i] || fXmetro1[i]>fMax[i] )
         in=0;
   }

   // calculate the log probabilities and compare old and new point
   double p0 = LogEvalSampling(fXmetro0); // old point
   double p1=0; // new point
   if(in)
      p1= LogEvalSampling(fXmetro1);

   // compare log probabilities
   int accept=0;
   if(in) {
      if(p1>=p0)
         accept=1;
      else {
         double r=log(fRandom->Rndm());
         if(r<p1-p0)
            accept=1;
      }
   }

   // fill the return point after the decision
   if(accept)
      for(int i=0;i<fNvar;i++) {
         fXmetro0[i]=fXmetro1[i];
         x[i]=fXmetro1[i];
      }
   else
      for(int i=0;i<fNvar;i++)
         x[i]=fXmetro0[i];

   fNmetro++;
}

void BCIntegrate::GetRandomVector ( std::vector< double > & x )

Fills a vector of random numbers between 0 and 1 into a vector

Parameters:

A

vector of doubles

Definition at line 1196 of file BCIntegrate.cxx.

{
   double * randx = new double[fNvar];

   fRandom->RndmArray(fNvar, randx);

   for(int i=0;i<fNvar;i++)
      x[i] = randx[i];

   delete[] randx;
   randx = 0;
}

void BCIntegrate::GetRandomVectorMetro ( std::vector< double > & x ) [virtual]

Definition at line 1210 of file BCIntegrate.cxx.

{
   double * randx = new double[fNvar];

   fRandom->RndmArray(fNvar, randx);

   for(int i=0;i<fNvar;i++)
      x[i] = 2.0 * (randx[i] - 0.5) * fMarkovChainStepSize;

   delete[] randx;
   randx = 0;
}

double BCIntegrate::GetRelativePrecision ( ) [inline]

Returns:: The relative precision for numerical integration

Definition at line 176 of file BCIntegrate.h.

         { return fRelativePrecision; };

BCIntegrate::BCSASchedule BCIntegrate::GetSASchedule ( ) [inline]

Returns:: The Simulated Annealing schedule

Definition at line 113 of file BCIntegrate.h.

         { return fSASchedule; };

double BCIntegrate::GetSAT0 ( ) [inline]

Returns the Simulated Annealing starting temperature.

Definition at line 218 of file BCIntegrate.h.

         { return fSAT0; }

double BCIntegrate::GetSATmin ( ) [inline]

Returns the Simulated Annealing threshhold temperature.

Definition at line 223 of file BCIntegrate.h.

         { return fSATmin; }

TTree* BCIntegrate::GetSATree ( ) [inline]

Definition at line 421 of file BCIntegrate.h.

         { return fTreeSA; };

void BCIntegrate::InitializeSATree ( )

Definition at line 1334 of file BCIntegrate.cxx.

{
   delete fTreeSA;
   fTreeSA = new TTree("SATree", "SATree");

   fTreeSA->Branch("Iteration",      &fSANIterations,   "iteration/I");
   fTreeSA->Branch("NParameters",    &fNvar,            "parameters/I");
   fTreeSA->Branch("Temperature",    &fSATemperature,   "temperature/D");
   fTreeSA->Branch("LogProbability", &fSALogProb,       "log(probability)/D");

   for (int i = 0; i < fNvar; ++i)
      fTreeSA->Branch(TString::Format("Parameter%i", i), &fSAx[i], TString::Format("parameter %i/D", i));
}

void BCIntegrate::InitMetro ( )

Initializes the Metropolis algorithm (for details see manual)

Definition at line 1071 of file BCIntegrate.cxx.

{
   fNmetro=0;

   if (fXmetro0.size() <= 0) {
      // start in the center of the phase space
      for(int i=0;i<fNvar;i++)
         fXmetro0.push_back((fMin[i]+fMax[i])/2.0);
   }

   fMarkovChainValue =  LogEval(fXmetro0);

   // run metropolis for a few times and dump the result... (to forget the initial position)
   std::vector <double> x;
   x.assign(fNvar,0.);

   for(int i=0;i<1000;i++)
      GetRandomPointMetro(x);

   fNmetro = 0;
}

double BCIntegrate::IntegralImportance ( std::vector< double > x )

Perfoms the importance sampling Monte Carlo integration. For details see documentation.

Parameters:

x

An initial point in parameter space

Returns:: The integral

Definition at line 549 of file BCIntegrate.cxx.

{
   // print debug information
   BCLog::OutDebug(Form("BCIntegrate::IntegralImportance. Integate over %i dimensions.", fNvar));

   // get total number of iterations
   double Niter = pow(fNiterPerDimension, fNvar);

   // print if more than 100,000 iterations
   if(Niter>1e5)
      BCLog::OutDebug(Form("BCIntegrate::IntegralImportance. Total number of iterations: %d.", Niter));

   // reset sum
   double sumI = 0;

   std::vector <double> randx;
   randx.assign(fNvar,0.);

   // prepare maximum value
   double pmax = 0.0;

   // loop over iterations
   for(int i = 0; i <= Niter; i++) {
      // get random point from sampling distribution
      GetRandomPointImportance(randx);

      // calculate probability at random point
      double val_f = Eval(randx);

      // calculate sampling distributions at that point
      double val_g = EvalSampling(randx);

      // add ratio to sum
      if (val_g > 0)
         sumI += val_f / val_g;

      // search for maximum probability
      if (val_f > pmax) {
         // set new maximum value
         pmax = val_f;

         // delete old best fit parameter values
         fBestFitParameters.clear();

         // write best fit parameters
         for(int i = 0; i < fNvar; i++)
            fBestFitParameters.push_back(randx.at(i));
      }

      // write intermediate results
      if((i+1)%100000 == 0)
         BCLog::OutDebug(Form("BCIntegrate::IntegralImportance. Iteration %d, integral: %d.", i+1, sumI/(i+1)));
   }

   // calculate integral
   double result=sumI/Niter;

   // print debug information
   BCLog::OutDebug(Form("BCIntegrate::IntegralImportance. Integral %f in %i iterations.", result, Niter));

   return result;
}

double BCIntegrate::IntegralMC	(	std::vector< double >	x,
		int *	varlist
	)

Perfoms the Monte Carlo integration. For details see documentation.

Parameters:

	x	An initial point in parameter space
	varlist	A list of variables

Returns:: The integral

Definition at line 382 of file BCIntegrate.cxx.

{
   SetVarList(varlist);
   return IntegralMC(x);
}

double BCIntegrate::IntegralMC ( std::vector< double > x )

Definition at line 389 of file BCIntegrate.cxx.

{
   // count the variables to integrate over
   int NvarNow=0;

   for(int i = 0; i < fNvar; i++)
      if(fVarlist[i])
         NvarNow++;

   // print to log
   BCLog::OutDebug(Form(" --> Running MC integation over %i dimensions.", NvarNow));
   BCLog::OutDebug(Form(" --> Maximum number of iterations: %i", GetNIterationsMax()));
   BCLog::OutDebug(Form(" --> Aimed relative precision:     %e", GetRelativePrecision()));

   // calculate D (the integrated volume)
   double D = 1.0;
   for(int i = 0; i < fNvar; i++)
      if (fVarlist[i])
         D = D * (fMax[i] - fMin[i]);

   // reset variables
   double pmax = 0.0;
   double sumW  = 0.0;
   double sumW2 = 0.0;
   double integral = 0.0;
   double variance = 0.0;
   double relprecision = 1.0;

   std::vector <double> randx;
   randx.assign(fNvar, 0.0);

   // reset number of iterations
   fNIterations = 0;

   // iterate while precision is not reached and the number of iterations is lower than maximum number of iterations
   while ((fRelativePrecision < relprecision && fNIterationsMax > fNIterations) || fNIterations < 10) {
      // increase number of iterations
      fNIterations++;

      // get random numbers
      GetRandomVector(randx);

      // scale random numbers
      for(int i = 0; i < fNvar; i++) {
         if(fVarlist[i])
            randx[i]=fMin[i]+randx[i]*(fMax[i]-fMin[i]);
         else
            randx[i]=x[i];
      }

      // evaluate function at sampled point
      double value = Eval(randx);

      // add value to sum and sum of squares
      sumW  += value;
      sumW2 += value * value;

      // search for maximum probability
      if (value > pmax) {
         // set new maximum value
         pmax = value;

         // delete old best fit parameter values
         fBestFitParameters.clear();

         // write best fit parameters
         for(int i = 0; i < fNvar; i++)
            fBestFitParameters.push_back(randx.at(i));
      }

      // calculate integral and variance
      integral = D * sumW / fNIterations;

      if (fNIterations%10000 == 0) {
         variance = (1.0 / double(fNIterations)) * (D * D * sumW2 / double(fNIterations) - integral * integral);
         double error = sqrt(variance);
         relprecision = error / integral;

         BCLog::Out(BCLog::debug, BCLog::debug,
            Form("BCIntegrate::IntegralMC. Iteration %i, integral: %e +- %e.", fNIterations, integral, error));
      }
   }

   fError = variance / fNIterations;

   // print to log
   BCLog::OutDebug(Form(" --> Result of integration:        %e +- %e   in %i iterations.", integral, sqrt(variance), fNIterations));
   BCLog::OutDebug(Form(" --> Obtained relative precision:  %e. ", sqrt(variance) / integral));

   return integral;
}

double BCIntegrate::IntegralMetro ( std::vector< double > x )

Perfoms the Metropolis Monte Carlo integration. For details see documentation.

Parameters:

x

An initial point in parameter space

Returns:: The integral

Definition at line 483 of file BCIntegrate.cxx.

{
   // print debug information
   BCLog::OutDebug(Form("BCIntegrate::IntegralMetro. Integate over %i dimensions.", fNvar));

   // get total number of iterations
   double Niter = pow(fNiterPerDimension, fNvar);

   // print if more than 100,000 iterations
   if(Niter>1e5)
      BCLog::OutDebug(Form(" --> Total number of iterations in Metropolis: %d.", Niter));

   // reset sum
   double sumI = 0;

   // prepare Metropolis algorithm
   std::vector <double> randx;
   randx.assign(fNvar,0.);
   InitMetro();

   // prepare maximum value
   double pmax = 0.0;

   // loop over iterations
   for(int i = 0; i <= Niter; i++) {
      // get random point from sampling distribution
      GetRandomPointSamplingMetro(randx);

      // calculate probability at random point
      double val_f = Eval(randx);

      // calculate sampling distributions at that point
      double val_g = EvalSampling(randx);

      // add ratio to sum
      if (val_g > 0)
         sumI += val_f / val_g;

      // search for maximum probability
      if (val_f > pmax) {
         // set new maximum value
         pmax = val_f;

         // delete old best fit parameter values
         fBestFitParameters.clear();

         // write best fit parameters
         for(int i = 0; i < fNvar; i++)
            fBestFitParameters.push_back(randx.at(i));
      }

      // write intermediate results
      if((i+1)%100000 == 0)
         BCLog::OutDebug(Form("BCIntegrate::IntegralMetro. Iteration %d, integral: %d.", i+1, sumI/(i+1)));
   }

   // calculate integral
   double result=sumI/Niter;

   // print debug information
   BCLog::OutDebug(Form(" --> Integral is %f in %i iterations.", result, Niter));

   return result;
}

double BCIntegrate::Integrate ( )

Does the integration over the un-normalized probability.

Returns:: The normalization value

Definition at line 350 of file BCIntegrate.cxx.

{
   std::vector <double> parameter;
   parameter.assign(fNvar, 0.0);

   switch(fIntegrationMethod)
   {
      case BCIntegrate::kIntMonteCarlo:
         return IntegralMC(parameter);

      case BCIntegrate::kIntMetropolis:
         return IntegralMetro(parameter);

      case BCIntegrate::kIntImportance:
         return IntegralImportance(parameter);

      case BCIntegrate::kIntCuba:
#ifdef HAVE_CUBA_H
         return CubaIntegrate();
#else
         BCLog::OutError("!!! This version of BAT is compiled without Cuba.");
         BCLog::OutError("    Use other integration methods or install Cuba and recompile BAT.");
#endif
   }

   BCLog::OutError(
         Form("BCIntegrate::Integrate : Invalid integration method: %d", fIntegrationMethod));

   return 0;
}

int BCIntegrate::IntegrateResetResults ( )

Reset all information on the best fit parameters.

Returns:: An error code

Definition at line 339 of file BCIntegrate.cxx.

{
   fBestFitParameters.clear();
   fBestFitParameterErrors.clear();
   fBestFitParametersMarginalized.clear();

   // no errors
   return 1;
}

double BCIntegrate::LogEval ( std::vector< double > x ) [virtual]

Evaluate the natural logarithm of the Eval function. For better numerical stability, this method should also be overloaded by the user.

Parameters:

x

The point in parameter space

Returns:: log(Eval(x))

Reimplemented from BCEngineMCMC.

Reimplemented in BCModel.

Definition at line 298 of file BCIntegrate.cxx.

{
   // this method should better also be overloaded
   return log(Eval(x));
}

double BCIntegrate::LogEvalSampling ( std::vector< double > x )

Evaluate the natural logarithm of the EvalSampling function. Method needs to be overloaded by the user.

Parameters:

x

The point in parameter space

Returns:: log(EvalSampling(x))

Definition at line 312 of file BCIntegrate.cxx.

{
   return log(EvalSampling(x));
}

TH1D * BCIntegrate::Marginalize ( BCParameter * parameter )

Performs the marginalization with respect to one parameter.

Parameters:

parameter

The parameter w.r.t. which the marginalization is performed

Returns:: A histogram which contains the marginalized probability distribution (normalized to 1)

Definition at line 613 of file BCIntegrate.cxx.

{
   BCLog::OutDebug(Form(" --> Marginalizing model wrt. parameter %s using method %d.", parameter->GetName().data(), fMarginalizationMethod));

   switch(fMarginalizationMethod)
   {
      case BCIntegrate::kMargMonteCarlo:
         return MarginalizeByIntegrate(parameter);

      case BCIntegrate::kMargMetropolis:
         return MarginalizeByMetro(parameter);
   }

   BCLog::OutWarning(Form("BCIntegrate::Marginalize. Invalid marginalization method: %d. Return 0.", fMarginalizationMethod));

   return 0;
}

TH2D * BCIntegrate::Marginalize	(	BCParameter *	parameter1,
		BCParameter *	parameter2
	)

Performs the marginalization with respect to two parameters.

Parameters:

	parameter1	The first parameter w.r.t. which the marginalization is performed
	parameter2	The second parameter w.r.t. which the marginalization is performed

Returns:: A histogram which contains the marginalized probability distribution (normalized to 1)

Definition at line 632 of file BCIntegrate.cxx.

{
   switch(fMarginalizationMethod)
   {
      case BCIntegrate::kMargMonteCarlo:
         return MarginalizeByIntegrate(parameter1,parameter2);

      case BCIntegrate::kMargMetropolis:
         return MarginalizeByMetro(parameter1,parameter2);
   }

   BCLog::OutWarning(Form("BCIntegrate::Marginalize. Invalid marginalization method: %d. Return 0.", fMarginalizationMethod));

   return 0;
}

int BCIntegrate::MarginalizeAllByMetro ( const char * name = "" )

Performs the marginalization with respect to every single parameter as well as with respect all combinations to two parameters using the Metropolis algorithm.

Parameters:

name

Basename for the histograms (e.g. model name)

Returns:: Total number of marginalized distributions

Definition at line 793 of file BCIntegrate.cxx.

{
   int niter=fNbins*fNbins*fNSamplesPer2DBin*fNvar;

   BCLog::OutDetail(Form(" --> Number of samples in Metropolis marginalization: %d.", niter));

   // define 1D histograms
   for(int i=0;i<fNvar;i++)
   {
      double hmin1 = fx->at(i)->GetLowerLimit();
      double hmax1 = fx->at(i)->GetUpperLimit();

      TH1D * h1 = new TH1D(
            TString::Format("h%s_%s", name, fx->at(i)->GetName().data()),"",
            fNbins, hmin1, hmax1);

      fHProb1D.push_back(h1);
   }

   // define 2D histograms
   for(int i=0;i<fNvar-1;i++)
      for(int j=i+1;j<fNvar;j++) {
         double hmin1 = fx->at(i)->GetLowerLimit();
         double hmax1 = fx->at(i)->GetUpperLimit();

         double hmin2 = fx->at(j)->GetLowerLimit();
         double hmax2 = fx->at(j)->GetUpperLimit();

         TH2D * h2 = new TH2D(
            TString::Format("h%s_%s_%s", name, fx->at(i)->GetName().data(), fx->at(j)->GetName().data()),"",
            fNbins, hmin1, hmax1,
            fNbins, hmin2, hmax2);

         fHProb2D.push_back(h2);
      }

   // get number of 2d distributions
   int nh2d = fHProb2D.size();

   BCLog::OutDetail(Form(" --> Marginalizing %d 1D distributions and %d 2D distributions.", fNvar, nh2d));

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

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

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

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

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

   // prepare Metro
   std::vector <double> randx;

   randx.assign(fNvar, 0.0);
   InitMetro();

   // reset counter
   fNacceptedMCMC = 0;

   // run Metro and fill histograms
   for(int i=0;i<=niter;i++) {
      GetRandomPointMetro(randx);

      // save this point to the markov chain in the ROOT file
      if (fFlagWriteMarkovChain) {
         fMCMCIteration = i;
         fMarkovChainTree->Fill();
      }

      for(int j=0;j<fNvar;j++)
         fHProb1D[j]->Fill( randx[j] );

      int ih=0;
      for(int j=0;j<fNvar-1;j++)
         for(int k=j+1;k<fNvar;k++) {
            fHProb2D[ih]->Fill(randx[j],randx[k]);
            ih++;
         }

      if((i+1)%100000==0)
         BCLog::OutDebug(Form("BCIntegrate::MarginalizeAllByMetro. %d samples done.",i+1));

      // function fitting

      if (fFitFunctionIndexX >= 0) {
         // loop over all possible x values ...
         if (fErrorBandContinuous) {
            double x = 0;

            for (int ix = 0; ix < fErrorBandNbinsX; ix++) {
               // calculate x
               x = fErrorBandXY->GetXaxis()->GetBinCenter(ix + 1);

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

               // fill histogram
               fErrorBandXY->Fill(x, y);
            }
         }

         // ... or evaluate at the data point x-values
         else {
            int ndatapoints = int(fErrorBandX.size());
            double x = 0;

            for (int ix = 0; ix < ndatapoints; ++ix) {
               // calculate x
               x = fErrorBandX.at(ix);

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

               // fill histogram
               fErrorBandXY->Fill(x, y);
            }
         }
      }
   }

   // normalize histograms
   for(int i=0;i<fNvar;i++)
      fHProb1D[i]->Scale( 1./fHProb1D[i]->Integral("width") );

   for (int i=0;i<nh2d;i++)
      fHProb2D[i]->Scale( 1./fHProb2D[i]->Integral("width") );

   if (fFitFunctionIndexX >= 0)
      fErrorBandXY->Scale(1.0/fErrorBandXY->Integral() * fErrorBandXY->GetNbinsX());

   if (fFitFunctionIndexX >= 0) {
      for (int ix = 1; ix <= fErrorBandNbinsX; ix++) {
         double sum = 0;

         for (int iy = 1; iy <= fErrorBandNbinsY; iy++)
            sum += fErrorBandXY->GetBinContent(ix, iy);

         for (int iy = 1; iy <= fErrorBandNbinsY; iy++) {
            double newvalue = fErrorBandXY->GetBinContent(ix, iy) / sum;
            fErrorBandXY->SetBinContent(ix, iy, newvalue);
         }
      }
   }

   BCLog::OutDebug(Form("BCIntegrate::MarginalizeAllByMetro done with %i trials and %i accepted trials. Efficiency is %f",fNmetro, fNacceptedMCMC, double(fNacceptedMCMC)/double(fNmetro)));

   return fNvar+nh2d;
}

TH1D * BCIntegrate::MarginalizeByIntegrate ( BCParameter * parameter )

Performs the marginalization with respect to one parameter using the simple Monte Carlo technique.

Parameters:

parameter

The parameter w.r.t. which the marginalization is performed

Returns:: A histogram which contains the marginalized probability distribution (normalized to 1)

Definition at line 649 of file BCIntegrate.cxx.

{
   // set parameter to marginalize
   ResetVarlist(1);
   int index = parameter->GetIndex();
   UnsetVar(index);

   // define histogram
   double hmin = parameter->GetLowerLimit();
   double hmax = parameter->GetUpperLimit();
   double hdx  = (hmax - hmin) / double(fNbins);
   TH1D * hist = new TH1D("hist","", fNbins, hmin, hmax);

   // integrate
   std::vector <double> randx;
   randx.assign(fNvar, 0.);

   for(int i=0;i<=fNbins;i++) {
      randx[index] = hmin + (double)i * hdx;

      double val = IntegralMC(randx);
      hist->Fill(randx[index], val);
   }

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

   return hist;
}

TH2D * BCIntegrate::MarginalizeByIntegrate	(	BCParameter *	parameter1,
		BCParameter *	parameter2
	)

Performs the marginalization with respect to two parameters using the simple Monte Carlo technique.

Parameters:

	parameter1	The first parameter w.r.t. which the marginalization is performed
	parameter2	The second parameter w.r.t. which the marginalization is performed

Returns:: A histogram which contains the marginalized probability distribution (normalized to 1)

Definition at line 680 of file BCIntegrate.cxx.

{
   // set parameter to marginalize
   ResetVarlist(1);
   int index1 = parameter1->GetIndex();
   UnsetVar(index1);
   int index2 = parameter2->GetIndex();
   UnsetVar(index2);

   // define histogram
   double hmin1 = parameter1->GetLowerLimit();
   double hmax1 = parameter1->GetUpperLimit();
   double hdx1  = (hmax1 - hmin1) / double(fNbins);

   double hmin2 = parameter2->GetLowerLimit();
   double hmax2 = parameter2->GetUpperLimit();
   double hdx2  = (hmax2 - hmin2) / double(fNbins);

   TH2D * hist = new TH2D(Form("hist_%s_%s", parameter1->GetName().data(), parameter2->GetName().data()),"",
         fNbins, hmin1, hmax1,
         fNbins, hmin2, hmax2);

   // integrate
   std::vector <double> randx;
   randx.assign(fNvar, 0.0);

   for(int i=0;i<=fNbins;i++) {
      randx[index1] = hmin1 + (double)i * hdx1;
      for(int j=0;j<=fNbins;j++) {
         randx[index2] = hmin2 + (double)j * hdx2;

         double val = IntegralMC(randx);
         hist->Fill(randx[index1],randx[index2], val);
      }
   }

   // normalize
   hist->Scale(1.0/hist->Integral("width"));

   return hist;
}

TH1D * BCIntegrate::MarginalizeByMetro ( BCParameter * parameter )

Performs the marginalization with respect to one parameter using the Metropolis algorithm.

Parameters:

parameter

The parameter w.r.t. which the marginalization is performed

Returns:: A histogram which contains the marginalized probability distribution (normalized to 1)

Definition at line 723 of file BCIntegrate.cxx.

{
   int niter = fMarkovChainNIterations;

   if (fMarkovChainAutoN == true)
      niter = fNbins*fNbins*fNSamplesPer2DBin*fNvar;

   BCLog::OutDetail(Form(" --> Number of samples in Metropolis marginalization: %d.", niter));

   // set parameter to marginalize
   int index = parameter->GetIndex();

   // define histogram
   double hmin = parameter->GetLowerLimit();
   double hmax = parameter->GetUpperLimit();
   TH1D * hist = new TH1D("hist","", fNbins, hmin, hmax);

   // prepare Metro
   std::vector <double> randx;
   randx.assign(fNvar, 0.0);
   InitMetro();

   for(int i=0;i<=niter;i++) {
      GetRandomPointMetro(randx);
      hist->Fill(randx[index]);
   }

   // normalize
   hist->Scale(1.0/hist->Integral("width"));

   return hist;
}

TH2D * BCIntegrate::MarginalizeByMetro	(	BCParameter *	parameter1,
		BCParameter *	parameter2
	)

Performs the marginalization with respect to two parameters using the Metropolis algorithm.

Parameters:

	parameter1	The first parameter w.r.t. which the marginalization is performed
	parameter2	The second parameter w.r.t. which the marginalization is performed

Returns:: A histogram which contains the marginalized probability distribution (normalized to 1)

Definition at line 757 of file BCIntegrate.cxx.

{
   int niter=fNbins*fNbins*fNSamplesPer2DBin*fNvar;

   // set parameter to marginalize
   int index1 = parameter1->GetIndex();
   int index2 = parameter2->GetIndex();

   // define histogram
   double hmin1 = parameter1->GetLowerLimit();
   double hmax1 = parameter1->GetUpperLimit();

   double hmin2 = parameter2->GetLowerLimit();
   double hmax2 = parameter2->GetUpperLimit();

   TH2D * hist = new TH2D(Form("hist_%s_%s", parameter1->GetName().data(), parameter2->GetName().data()),"",
         fNbins, hmin1, hmax1,
         fNbins, hmin2, hmax2);

   // prepare Metro
   std::vector <double> randx;
   randx.assign(fNvar, 0.0);
   InitMetro();

   for(int i=0;i<=niter;i++) {
      GetRandomPointMetro(randx);
      hist->Fill(randx[index1],randx[index2]);
   }

   // normalize
   hist->Scale(1.0/hist->Integral("width"));

   return hist;
}

void BCIntegrate::MCMCFillErrorBand ( ) [private]

Definition at line 1961 of file BCIntegrate.cxx.

{
   if (!fFillErrorBand)
      return;

   // function fitting
   if (fFitFunctionIndexX < 0)
      return;

   // loop over all possible x values ...
   if (fErrorBandContinuous) {
      double x = 0;
      for (int ix = 0; ix < fErrorBandNbinsX; ix++) {
         // calculate x
         x = fErrorBandXY->GetXaxis()->GetBinCenter(ix + 1);

         // calculate y
         std::vector <double> xvec;
         xvec.push_back(x);

         // loop over all chains
         for (int ichain = 0; ichain < MCMCGetNChains(); ++ichain) {
            // calculate y
            double y = FitFunction(xvec, MCMCGetx(ichain));

            // fill histogram
            fErrorBandXY->Fill(x, y);
         }

         xvec.clear();
      }
   }
   // ... or evaluate at the data point x-values
   else {
      int ndatapoints = int(fErrorBandX.size());
      double x = 0;

      for (int ix = 0; ix < ndatapoints; ++ix) {
         // calculate x
         x = fErrorBandX.at(ix);

         // calculate y
         std::vector <double> xvec;
         xvec.push_back(x);

         // loop over all chains
         for (int ichain = 0; ichain < MCMCGetNChains(); ++ichain) {
            // calculate y
            double y = FitFunction(xvec, MCMCGetx(ichain));

            // fill histogram
            fErrorBandXY->Fill(x, y);
         }

         xvec.clear();
      }
   }
}

void BCIntegrate::MCMCIterationInterface ( ) [private, virtual]

Reimplemented from BCEngineMCMC.

Definition at line 1948 of file BCIntegrate.cxx.

{
   // what's within this method will be executed
   // for every iteration of the MCMC

   // fill error band
   MCMCFillErrorBand();

   // do user defined stuff
   MCMCUserIterationInterface();
}

virtual void BCIntegrate::MCMCUserIterationInterface ( ) [inline, virtual]

Reimplemented in BCGoFTest, and BCTemplateFitter.

Definition at line 717 of file BCIntegrate.h.

{};

void BCIntegrate::ResetVarlist ( int v )

Sets all values of the variable list to a particular value The value

Definition at line 284 of file BCIntegrate.cxx.

{
   for(int i=0;i<fNvar;i++)
      fVarlist[i]=v;
}

double BCIntegrate::SAHelperGetRadialCauchy ( )

Generates the radial part of a n-dimensional Cauchy distribution. Helper function for Cauchy annealing.

Definition at line 1632 of file BCIntegrate.cxx.

{
   // theta is sampled from a rather complicated distribution,
   // so first we create a lookup table with 10000 random numbers
   // once and then, each time we need a new random number,
   // we just look it up in the table.
   double theta;

   // static vectors for theta-sampling-map
   static double map_u[10001];
   static double map_theta[10001];
   static bool initialized = false;
   static int map_dimension = 0;

   // is the lookup-table already initialized? if not, do it!
   if (!initialized || map_dimension != fNvar) {
      double init_theta;
      double beta = SAHelperSinusToNIntegral(fNvar - 1, 1.57079632679);

      for (int i = 0; i <= 10000; i++) {
         init_theta = 3.14159265 * (double)i / 5000.;
         map_theta[i] = init_theta;

         map_u[i] = SAHelperSinusToNIntegral(fNvar - 1, init_theta) / beta;
      }

      map_dimension = fNvar;
      initialized = true;
   } // initializing is done.

   // generate uniform random number for sampling
   double u = fRandom->Rndm();

   // Find the two elements just greater than and less than u
   // using a binary search (O(log(N))).
   int lo = 0;
   int up = 10000;
   int mid;

   while (up != lo) {
      mid = ((up - lo + 1) / 2) + lo;

      if (u >= map_u[mid])
         lo = mid;
      else
         up = mid - 1;
   }
   up++;

   // perform linear interpolation:
   theta = map_theta[lo] + (u - map_u[lo]) / (map_u[up] - map_u[lo]) * (map_theta[up] - map_theta[lo]);

   return tan(theta);
}

std::vector< double > BCIntegrate::SAHelperGetRandomPointOnHypersphere ( )

Generates a uniform distributed random point on the surface of a fNvar-dimensional Hypersphere. Used as a helper to generate proposal points for Cauchy annealing.

Definition at line 1587 of file BCIntegrate.cxx.

{
   std::vector<double> rand_point(fNvar);

   // This method can only be called with fNvar >= 2 since the 1-dim case
   // is already hard wired into the Cauchy annealing proposal function.
   // To speed things up, hard-code fast method for 2 and dimensions.
   // The algorithm for 2D can be found at
   // http://mathworld.wolfram.com/CirclePointPicking.html
   // For 3D just using ROOT's algorithm.
   if (fNvar == 2) {
      double x1, x2, s;
      do {
         x1 = fRandom->Rndm() * 2. - 1.;
         x2 = fRandom->Rndm() * 2. - 1.;
         s = x1*x1 + x2*x2;
      }
      while (s >= 1);

      rand_point[0] = (x1*x1 - x2*x2) / s;
      rand_point[1] = (2.*x1*x2) / s;
   }
   else if (fNvar == 3) {
      fRandom->Sphere(rand_point[0], rand_point[1], rand_point[2], 1.0);
   }
   else {
      double s = 0.,
         gauss_num;

      for (int i = 0; i < fNvar; i++) {
         gauss_num = fRandom->Gaus();
         rand_point[i] = gauss_num;
         s += gauss_num * gauss_num;
      }
      s = sqrt(s);

      for (int i = 0; i < fNvar; i++)
         rand_point[i] = rand_point[i] / s;
   }

   return rand_point;
}

double BCIntegrate::SAHelperSinusToNIntegral	(	int	dim,
		double	theta
	)

Returns the Integral of sin^dim from 0 to theta. Helper function needed for generating Cauchy distributions.

Definition at line 1689 of file BCIntegrate.cxx.

{
   if (dim < 1)
      return theta;
   else if (dim == 1)
      return (1. - cos(theta));
   else if (dim == 2)
      return 0.5 * (theta - sin(theta) * cos(theta));
   else if (dim == 3)
      return (2. - sin(theta) * sin(theta) * cos(theta) - 2. * cos(theta)) / 3.;
   else
      return - pow(sin(theta), (double)(dim - 1)) * cos(theta) / (double)dim
         + (double)(dim - 1) / (double)dim
         * SAHelperSinusToNIntegral(dim - 2, theta);
}

void BCIntegrate::SAInitialize ( )

Definition at line 2021 of file BCIntegrate.cxx.

{
   fSAx.clear();
   fSAx.assign(fNvar, 0.0);
}

double BCIntegrate::SATemperature ( double t )

Temperature annealing schedule for use with Simulated Annealing. Delegates to the appropriate method according to fSASchedule.

Parameters:

t

iterator for lowering the temperature over time.

Definition at line 1479 of file BCIntegrate.cxx.

{
   // do we have Cauchy (default), Boltzmann or custom annealing schedule?
   if (fSASchedule == BCIntegrate::kSABoltzmann)
      return SATemperatureBoltzmann(t);
   else if (fSASchedule == BCIntegrate::kSACauchy)
      return SATemperatureCauchy(t);
   else
      return SATemperatureCustom(t);
}

double BCIntegrate::SATemperatureBoltzmann ( double t )

Temperature annealing schedule for use with Simulated Annealing. This method is used for Boltzmann annealing schedule.

Parameters:

t

iterator for lowering the temperature over time.

Definition at line 1492 of file BCIntegrate.cxx.

{
   return fSAT0 / log((double)(t + 1));
}

double BCIntegrate::SATemperatureCauchy ( double t )

Temperature annealing schedule for use with Simulated Annealing. This method is used for Cauchy annealing schedule.

Parameters:

t

iterator for lowering the temperature over time.

Definition at line 1499 of file BCIntegrate.cxx.

{
   return fSAT0 / (double)t;
}

double BCIntegrate::SATemperatureCustom ( double t ) [virtual]

Temperature annealing schedule for use with Simulated Annealing. This is a virtual method to be overridden by a user-defined custom temperature schedule.

Parameters:

t

iterator for lowering the temperature over time.

Definition at line 1506 of file BCIntegrate.cxx.

{
   BCLog::OutError("BCIntegrate::SATemperatureCustom : No custom temperature schedule defined");
   return 0.;
}

void BCIntegrate::SetDataPointLowerBoundaries ( BCDataPoint * datasetlowerboundaries ) [inline]

Sets the data point containing the lower boundaries of possible data values

Definition at line 338 of file BCIntegrate.h.

         { fDataPointLowerBoundaries = datasetlowerboundaries; };

void BCIntegrate::SetDataPointLowerBoundary	(	int	index,
		double	lowerboundary
	)			`[inline]`

Sets the lower boundary of possible data values for a particular variable

Definition at line 350 of file BCIntegrate.h.

         { fDataPointLowerBoundaries -> SetValue(index, lowerboundary); };

void BCIntegrate::SetDataPointUpperBoundaries ( BCDataPoint * datasetupperboundaries ) [inline]

Sets the data point containing the upper boundaries of possible data values

Definition at line 344 of file BCIntegrate.h.

         { fDataPointUpperBoundaries = datasetupperboundaries; };

void BCIntegrate::SetDataPointUpperBoundary	(	int	index,
		double	upperboundary
	)			`[inline]`

Sets the upper boundary of possible data values for a particular variable

Definition at line 356 of file BCIntegrate.h.

         { fDataPointUpperBoundaries -> SetValue(index, upperboundary); };

void BCIntegrate::SetErrorBandHisto ( TH2D * h ) [inline]

Sets errorband histogram

Definition at line 398 of file BCIntegrate.h.

         { fErrorBandXY = h; };

void BCIntegrate::SetFillErrorBand ( bool flag = true ) [inline]

Definition at line 310 of file BCIntegrate.h.

         { fFillErrorBand=flag; };

void BCIntegrate::SetFitFunctionIndexX ( int index ) [inline]

Sets index of the x values in function fits.

Parameters:

index

Index of the x values

Definition at line 322 of file BCIntegrate.h.

         { fFitFunctionIndexX = index; };

void BCIntegrate::SetFitFunctionIndexY ( int index ) [inline]

Sets index of the y values in function fits.

Parameters:

index

Index of the y values

Definition at line 328 of file BCIntegrate.h.

         { fFitFunctionIndexY = index; };

void BCIntegrate::SetFitFunctionIndices	(	int	indexx,
		int	indexy
	)			`[inline]`

Definition at line 331 of file BCIntegrate.h.

         { this -> SetFitFunctionIndexX(indexx);
            this -> SetFitFunctionIndexY(indexy); };

void BCIntegrate::SetFlagIgnorePrevOptimization ( bool flag ) [inline]

Definition at line 235 of file BCIntegrate.h.

         { fFlagIgnorePrevOptimization = flag; };

void BCIntegrate::SetFlagWriteSAToFile ( bool flag ) [inline]

Definition at line 411 of file BCIntegrate.h.

         { fFlagWriteSAToFile = flag; };

void BCIntegrate::SetIntegrationMethod ( BCIntegrate::BCIntegrationMethod method )

Parameters:

method

The integration method

Definition at line 327 of file BCIntegrate.cxx.

{
#ifdef HAVE_CUBA_H
   fIntegrationMethod = method;
#else
   BCLog::OutWarning("!!! This version of BAT is compiled without Cuba.");
   BCLog::OutWarning("    Monte Carlo Sampled Mean integration method will be used.");
   BCLog::OutWarning("    To be able to use Cuba integration, install Cuba and recompile BAT.");
#endif
}

void BCIntegrate::SetMarginalizationMethod ( BCIntegrate::BCMarginalizationMethod method ) [inline]

Parameters:

method

The marginalization method

Definition at line 257 of file BCIntegrate.h.

         { fMarginalizationMethod = method; };

void BCIntegrate::SetMarkovChainAutoN ( bool flag ) [inline]

Sets a flag for automatically calculating the number of iterations

Definition at line 389 of file BCIntegrate.h.

         { fMarkovChainAutoN = flag; };

void BCIntegrate::SetMarkovChainInitialPosition ( std::vector< double > position )

Sets the initial position for the Markov chain

Definition at line 187 of file BCIntegrate.cxx.

{
   int n = position.size();

   fXmetro0.clear();

   for (int i = 0; i < n; ++i)
      fXmetro0.push_back(position.at(i));
}

void BCIntegrate::SetMarkovChainNIterations ( int niterations ) [inline]

Sets the number of iterations in the markov chain

Definition at line 383 of file BCIntegrate.h.

         { fMarkovChainNIterations = niterations;
           fMarkovChainAutoN = false; }

void BCIntegrate::SetMarkovChainStepSize ( double stepsize ) [inline]

Sets the step size for Markov chains

Definition at line 378 of file BCIntegrate.h.

         { fMarkovChainStepSize = stepsize; };

void BCIntegrate::SetMarkovChainTree ( TTree * tree ) [inline]

Sets the ROOT tree containing the Markov chain

Definition at line 369 of file BCIntegrate.h.

         { fMarkovChainTree = tree; };

void BCIntegrate::SetMinuitArlist ( double * arglist ) [inline]

Definition at line 231 of file BCIntegrate.h.

         { fMinuitArglist[0] = arglist[0];
           fMinuitArglist[1] = arglist[1]; };

void BCIntegrate::SetMode ( std::vector< double > mode )

Sets mode

Definition at line 1707 of file BCIntegrate.cxx.

{
   if((int)mode.size() == fNvar) {
      fBestFitParameters.clear();
      for (int i = 0; i < fNvar; i++)
         fBestFitParameters.push_back(mode[i]);
   }
}

void BCIntegrate::SetNbins	(	int	nbins,
		int	index = `-1`
	)

Parameters:

	n	Number of bins per dimension for the marginalized distributions. Default is 100. Minimum number allowed is 2. Set the number of bins for the marginalized distribution of a parameter.
	nbins	Number of bins (default = 100)
	index	Index of the parameter.

Definition at line 216 of file BCIntegrate.cxx.

{
   if (fNvar == 0)
      return;

   // check if index is in range
   if (index >= fNvar) {
      BCLog::OutWarning("BCIntegrate::SetNbins : Index out of range.");
      return;
   }
   // set for all parameters at once
   else if (index < 0) {
      for (int i = 0; i < fNvar; ++i)
         SetNbins(nbins, i);
      return;
   }

   // sanity check for nbins
   if (nbins <= 0)
      nbins = 10;

   fMCMCH1NBins[index] = nbins;

   return;

//    if(n<2) {
//       BCLog::OutWarning("BCIntegrate::SetNbins. Number of bins less than 2 makes no sense. Setting to 2.");
//       n=2;
//    }
//    MCMCSetH1NBins(n, -1);

   //   fNbins=n;

   //   fMCMCH1NBins = n;
   //   fMCMCH2NBinsX = n;
   //   fMCMCH2NBinsY = n;
}

void BCIntegrate::SetNIterationsMax ( int niterations ) [inline]

Parameters:

niterations

The maximum number of iterations for Monte Carlo integration

Definition at line 283 of file BCIntegrate.h.

         { fNIterationsMax = niterations; };

void BCIntegrate::SetNiterationsPerDimension ( int niterations ) [inline]

Parameters:

niterations

Number of iterations per dimension for Monte Carlo integration.

Definition at line 272 of file BCIntegrate.h.

         { fNiterPerDimension = niterations; };

void BCIntegrate::SetNSamplesPer2DBin ( int n ) [inline]

Parameters:

n

Number of samples per 2D bin per variable in the Metropolis marginalization. Default is 100.

Definition at line 278 of file BCIntegrate.h.

         { fNSamplesPer2DBin = n; };

void BCIntegrate::SetOptimizationMethod ( BCIntegrate::BCOptimizationMethod method ) [inline]

Parameters:

method

The mode finding method

Definition at line 262 of file BCIntegrate.h.

         { fOptimizationMethod = method; };

void BCIntegrate::SetParameters ( BCParameterSet * par )

Parameters:

par

The parameter set which gets translated into array needed for the Monte Carlo integration

Definition at line 147 of file BCIntegrate.cxx.

{
   DeleteVarList();

   fx = par;
   fNvar = fx->size();
   fMin = new double[fNvar];
   fMax = new double[fNvar];
   fVarlist = new int[fNvar];

   ResetVarlist(1);

   for(int i=0;i<fNvar;i++) {
      fMin[i]=(fx->at(i))->GetLowerLimit();
      fMax[i]=(fx->at(i))->GetUpperLimit();
   }

   fXmetro0.clear();
   for(int i=0;i<fNvar;i++)
      fXmetro0.push_back((fMin[i]+fMax[i])/2.0);

   fXmetro1.clear();
   fXmetro1.assign(fNvar,0.);

   fMCMCBoundaryMin.clear();
   fMCMCBoundaryMax.clear();

   for(int i=0;i<fNvar;i++) {
      fMCMCBoundaryMin.push_back(fMin[i]);
      fMCMCBoundaryMax.push_back(fMax[i]);
      fMCMCFlagsFillHistograms.push_back(true);
   }

   for (int i = int(fMCMCH1NBins.size()); i<fNvar; ++i)
      fMCMCH1NBins.push_back(100);

   fMCMCNParameters = fNvar;
}

void BCIntegrate::SetRelativePrecision ( double relprecision ) [inline]

Parameters:

relprecision

The relative precision envisioned for Monte Carlo integration

Definition at line 289 of file BCIntegrate.h.

         { fRelativePrecision = relprecision; };

void BCIntegrate::SetSASchedule ( BCIntegrate::BCSASchedule schedule ) [inline]

Parameters:

method

The Simulated Annealing schedule

Definition at line 267 of file BCIntegrate.h.

         { fSASchedule = schedule; };

void BCIntegrate::SetSAT0 ( double T0 ) [inline]

Parameters:

T0

new value for Simulated Annealing starting temperature.

Definition at line 403 of file BCIntegrate.h.

         { fSAT0 = T0; }

void BCIntegrate::SetSATmin ( double Tmin ) [inline]

Parameters:

Tmin

new value for Simulated Annealing threshold temperature.

Definition at line 408 of file BCIntegrate.h.

         { fSATmin = Tmin; }

void BCIntegrate::SetSATree ( TTree * tree ) [inline]

Definition at line 416 of file BCIntegrate.h.

         { fTreeSA = tree; };

void BCIntegrate::SetVar ( int index ) [inline]

Parameters:

index

The index of the variable to be set

Definition at line 249 of file BCIntegrate.h.

{fVarlist[index]=1;};

void BCIntegrate::SetVarList ( int * varlist )

Parameters:

varlist

A list of parameters

Definition at line 277 of file BCIntegrate.cxx.

{
   for(int i=0;i<fNvar;i++)
      fVarlist[i]=varlist[i];
}

void BCIntegrate::UnsetFillErrorBand ( ) [inline]

Definition at line 316 of file BCIntegrate.h.

         { fFillErrorBand=false; };

void BCIntegrate::UnsetVar ( int index ) [inline]

Set value of a particular integration variable to 0.

Parameters:

index

The index of the variable

Definition at line 445 of file BCIntegrate.h.

         { fVarlist[index] = 0; };

void BCIntegrate::WriteMarkovChain ( bool flag ) [inline]

Flag for writing Markov chain to ROOT file (true) or not (false)

Definition at line 362 of file BCIntegrate.h.

         { fFlagWriteMarkovChain = flag;
           fMCMCFlagWriteChainToFile = flag;
           fMCMCFlagWritePreRunToFile = flag; };

Member Data Documentation

std::vector<double> BCIntegrate::fBestFitParameterErrors [protected]

Definition at line 848 of file BCIntegrate.h.

std::vector<double> BCIntegrate::fBestFitParameters [protected]

A vector of best fit parameters estimated from the global probability and the estimate on their uncertainties

Definition at line 847 of file BCIntegrate.h.

std::vector<double> BCIntegrate::fBestFitParametersMarginalized [protected]

A vector of best fit parameters estimated from the marginalized probability

Definition at line 852 of file BCIntegrate.h.

std::vector<bool> BCIntegrate::fDataFixedValues [protected]

Definition at line 842 of file BCIntegrate.h.

BCDataPoint* BCIntegrate::fDataPointLowerBoundaries [protected]

data point containing the lower boundaries of possible data values

Definition at line 836 of file BCIntegrate.h.

BCDataPoint* BCIntegrate::fDataPointUpperBoundaries [protected]

data point containing the upper boundaries of possible data values

Definition at line 840 of file BCIntegrate.h.

double BCIntegrate::fError [private]

The uncertainty in the most recent Monte Carlo integration

Definition at line 784 of file BCIntegrate.h.

bool BCIntegrate::fErrorBandContinuous [protected]

A flag for single point evaluation of the error "band"

Definition at line 873 of file BCIntegrate.h.

int BCIntegrate::fErrorBandNbinsX [protected]

Number of X bins of the error band histogram

Definition at line 882 of file BCIntegrate.h.

int BCIntegrate::fErrorBandNbinsY [protected]

Nnumber of Y bins of the error band histogram

Definition at line 886 of file BCIntegrate.h.

std::vector<double> BCIntegrate::fErrorBandX [protected]

Definition at line 874 of file BCIntegrate.h.

TH2D* BCIntegrate::fErrorBandXY [protected]

The error band histogram

Definition at line 878 of file BCIntegrate.h.

bool BCIntegrate::fFillErrorBand [protected]

Flag whether or not to fill the error band

Definition at line 864 of file BCIntegrate.h.

int BCIntegrate::fFitFunctionIndexX [protected]

The indeces for function fits

Definition at line 868 of file BCIntegrate.h.

int BCIntegrate::fFitFunctionIndexY [protected]

Definition at line 869 of file BCIntegrate.h.

double BCIntegrate::fFlagIgnorePrevOptimization [protected]

Flag for ignoring older results of minimization

Definition at line 897 of file BCIntegrate.h.

bool BCIntegrate::fFlagWriteMarkovChain [protected]

Flag for writing Markov chain to file

Definition at line 901 of file BCIntegrate.h.

bool BCIntegrate::fFlagWriteSAToFile [protected]

Flag deciding whether to write SA to file or not.

Definition at line 925 of file BCIntegrate.h.

std::vector<TH1D *> BCIntegrate::fHProb1D [protected]

Vector of TH1D histograms for marginalized probability distributions

Definition at line 856 of file BCIntegrate.h.

std::vector<TH2D *> BCIntegrate::fHProb2D [protected]

Vector of TH2D histograms for marginalized probability distributions

Definition at line 860 of file BCIntegrate.h.

BCIntegrate::BCIntegrationMethod BCIntegrate::fIntegrationMethod [private]

Current integration method

Definition at line 751 of file BCIntegrate.h.

BCIntegrate::BCMarginalizationMethod BCIntegrate::fMarginalizationMethod [private]

Current marginalization method

Definition at line 755 of file BCIntegrate.h.

bool BCIntegrate::fMarkovChainAutoN [protected]

Definition at line 832 of file BCIntegrate.h.

int BCIntegrate::fMarkovChainNIterations [protected]

Definition at line 830 of file BCIntegrate.h.

double BCIntegrate::fMarkovChainStepSize [protected]

Step size in the Markov chain relative to min and max

Definition at line 828 of file BCIntegrate.h.

TTree* BCIntegrate::fMarkovChainTree [protected]

ROOT tree containing the Markov chain

Definition at line 905 of file BCIntegrate.h.

double BCIntegrate::fMarkovChainValue [private]

A double containing the log likelihood value at the point fXmetro1

Definition at line 801 of file BCIntegrate.h.

double* BCIntegrate::fMax [private]

Array containing the upper boundaries of the variables to integrate over.

Definition at line 739 of file BCIntegrate.h.

int BCIntegrate::fMCMCIteration [protected]

Iteration of the MCMC

Definition at line 909 of file BCIntegrate.h.

double* BCIntegrate::fMin [private]

Array containing the lower boundaries of the variables to integrate over.

Definition at line 735 of file BCIntegrate.h.

TMinuit* BCIntegrate::fMinuit [protected]

Minuit

Definition at line 890 of file BCIntegrate.h.

double BCIntegrate::fMinuitArglist[2] [protected]

Definition at line 892 of file BCIntegrate.h.

int BCIntegrate::fMinuitErrorFlag [protected]

Definition at line 893 of file BCIntegrate.h.

int BCIntegrate::fNacceptedMCMC [private]

Definition at line 789 of file BCIntegrate.h.

int BCIntegrate::fNbins [protected]

Number of bins per dimension for the marginalized distributions

Definition at line 819 of file BCIntegrate.h.

int BCIntegrate::fNIterations [private]

Number of iterations in the most recent Monte Carlo integation

Definition at line 776 of file BCIntegrate.h.

int BCIntegrate::fNIterationsMax [private]

Maximum number of iterations

Definition at line 772 of file BCIntegrate.h.

int BCIntegrate::fNiterPerDimension [private]

Number of iteration per dimension for Monte Carlo integration.

Definition at line 747 of file BCIntegrate.h.

int BCIntegrate::fNmetro [private]

The number of iterations in the Metropolis integration

Definition at line 788 of file BCIntegrate.h.

int BCIntegrate::fNSamplesPer2DBin [protected]

Number of samples per 2D bin per variable in the Metropolis marginalization.

Definition at line 824 of file BCIntegrate.h.

int BCIntegrate::fNvar [protected]

Number of variables to integrate over.

Definition at line 815 of file BCIntegrate.h.

BCIntegrate::BCOptimizationMethod BCIntegrate::fOptimizationMethod [private]

Current mode finding method

Definition at line 759 of file BCIntegrate.h.

BCIntegrate::BCOptimizationMethod BCIntegrate::fOptimizationMethodMode [private]

Method with which the global mode was found (can differ from fOptimization method in case more than one algorithm is used).

Definition at line 764 of file BCIntegrate.h.

double BCIntegrate::fRelativePrecision [private]

Relative precision aimed at in the Monte Carlo integation

Definition at line 780 of file BCIntegrate.h.

double BCIntegrate::fSALogProb [protected]

Definition at line 929 of file BCIntegrate.h.

int BCIntegrate::fSANIterations [protected]

Definition at line 927 of file BCIntegrate.h.

BCIntegrate::BCSASchedule BCIntegrate::fSASchedule [private]

Current Simulated Annealing schedule

Definition at line 768 of file BCIntegrate.h.

double BCIntegrate::fSAT0 [protected]

Starting temperature for Simulated Annealing

Definition at line 913 of file BCIntegrate.h.

double BCIntegrate::fSATemperature [protected]

Definition at line 928 of file BCIntegrate.h.

double BCIntegrate::fSATmin [protected]

Minimal/Threshold temperature for Simulated Annealing

Definition at line 917 of file BCIntegrate.h.

std::vector<double> BCIntegrate::fSAx [protected]

Definition at line 930 of file BCIntegrate.h.

TTree* BCIntegrate::fTreeSA [protected]

Tree for the Simulated Annealing

Definition at line 921 of file BCIntegrate.h.

int* BCIntegrate::fVarlist [private]

List of variables containing a flag whether to integrate over them or not.

Definition at line 743 of file BCIntegrate.h.

BCParameterSet* BCIntegrate::fx [private]

Set of parameters for the integration.

Definition at line 731 of file BCIntegrate.h.

std::vector<double> BCIntegrate::fXmetro0 [private]

A vector of points in parameter space used for the Metropolis algorithm

Definition at line 793 of file BCIntegrate.h.

std::vector<double> BCIntegrate::fXmetro1 [private]

A vector of points in parameter space used for the Metropolis algorithm

Definition at line 797 of file BCIntegrate.h.

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

BCIntegrate Class Reference

Public Types

Public Member Functions

Protected Attributes

Private Member Functions

Private Attributes

Member functions (miscellaneous methods)

Detailed Description

Member Enumeration Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation