BCIntegrate Class Reference

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

#include <BCIntegrate.h>

Inherits BCEngineMCMC.

Inherited by BCModel.

Collaboration diagram for BCIntegrate:

[legend]

List of all members.

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 (BCParameterSet *par)
	BCIntegrate ()
virtual	~BCIntegrate ()
Member functions (get)
double	GetError ()
BCIntegrate::BCIntegrationMethod	GetIntegrationMethod ()
BCIntegrate::BCMarginalizationMethod	GetMarginalizationMethod ()
std::vector< double > *	GetMarkovChainPoint ()
TTree *	GetMarkovChainTree ()
double *	GetMarkovChainValue ()
int *	GetMCMCIteration ()
TMinuit *	GetMinuit ()
int	GetMinuitErrorFlag ()
int	GetNbins ()
int	GetNIterations ()
int	GetNIterationsMax ()
int	GetNiterationsPerDimension ()
int	GetNSamplesPer2DBin ()
int	GetNvar ()
BCIntegrate::BCOptimizationMethod	GetOptimizationMethod ()
BCIntegrate::BCOptimizationMethod	GetOptimizationMethodMode ()
double	GetRandomPoint (std::vector< double > &x)
double	GetRandomPointImportance (std::vector< double > &x)
void	GetRandomPointMetro (std::vector< double > &x)
void	GetRandomPointSamplingMetro (std::vector< double > &x)
void	GetRandomVector (std::vector< double > &x)
virtual void	GetRandomVectorMetro (std::vector< double > &x)
double	GetRelativePrecision ()
double *	GetSAP2LogProb ()
int *	GetSAP2NIterations ()
double *	GetSAP2Temperature ()
std::vector< double > *	GetSAP2x ()
BCIntegrate::BCSASchedule	GetSASchedule ()
double	GetSAT0 ()
double	GetSATmin ()
Member functions (set)
void	SetDataPointLowerBoundaries (BCDataPoint *datasetlowerboundaries)
void	SetDataPointLowerBoundary (int index, double lowerboundary)
void	SetDataPointUpperBoundaries (BCDataPoint *datasetupperboundaries)
void	SetDataPointUpperBoundary (int index, double upperboundary)
void	SetErrorBandHisto (TH2D *h)
void	SetFillErrorBand (bool flag=true)
void	SetFitFunctionIndexX (int index)
void	SetFitFunctionIndexY (int index)
void	SetFitFunctionIndices (int indexx, int indexy)
void	SetFlagIgnorePrevOptimization (bool flag)
void	SetFlagWriteSAToFile (bool flag)
void	SetIntegrationMethod (BCIntegrate::BCIntegrationMethod method)
void	SetMarginalizationMethod (BCIntegrate::BCMarginalizationMethod method)
void	SetMarkovChainAutoN (bool flag)
void	SetMarkovChainInitialPosition (std::vector< double > position)
void	SetMarkovChainNIterations (int niterations)
void	SetMarkovChainStepSize (double stepsize)
void	SetMarkovChainTree (TTree *tree)
void	SetMinuitArlist (double *arglist)
void	SetMode (std::vector< double > mode)
void	SetNbins (int nbins, int index=-1)
void	SetNIterationsMax (int niterations)
void	SetNiterationsPerDimension (int niterations)
void	SetNSamplesPer2DBin (int n)
void	SetOptimizationMethod (BCIntegrate::BCOptimizationMethod method)
void	SetParameters (BCParameterSet *par)
void	SetRelativePrecision (double relprecision)
void	SetSASchedule (BCIntegrate::BCSASchedule schedule)
void	SetSAT0 (double T0)
void	SetSATmin (double Tmin)
void	SetSATree (TTree *tree)
void	SetVar (int index)
void	SetVarList (int *varlist)
void	UnsetFillErrorBand ()
void	WriteMarkovChain (bool flag)
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
TRandom3 *	fRandom
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)
double	CubaIntegrate ()
double	CubaIntegrate (int method, std::vector< double > parameters_double, std::vector< double > parameters_int)
void	DeleteVarList ()
virtual double	Eval (std::vector< double > x)
double	EvalPrint (std::vector< double > x)
virtual double	EvalSampling (std::vector< double > x)
void	FindModeMCMC ()
virtual void	FindModeMinuit (std::vector< double > start=std::vector< double >(0), int printlevel=1)
void	FindModeSA (std::vector< double > start=std::vector< double >(0))
virtual double	FitFunction (std::vector< double > x, std::vector< double > parameters)
TH1D *	GetH1D (int parIndex)
TH2D *	GetH2D (int parIndex1, int parIndex2)
int	GetH2DIndex (int parIndex1, int parIndex2)
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)
void	InitMetro ()
double	IntegralImportance (std::vector< double > x)
double	IntegralMC (std::vector< double > x)
double	IntegralMC (std::vector< double > x, int *varlist)
double	IntegralMetro (std::vector< double > x)
double	Integrate ()
virtual double	LogEval (std::vector< double > x)
double	LogEvalSampling (std::vector< double > x)
TH2D *	Marginalize (BCParameter *parameter1, BCParameter *parameter2)
TH1D *	Marginalize (BCParameter *parameter)
int	MarginalizeAllByMetro (const char *name)
TH2D *	MarginalizeByIntegrate (BCParameter *parameter1, BCParameter *parameter2)
TH1D *	MarginalizeByIntegrate (BCParameter *parameter)
TH2D *	MarginalizeByMetro (BCParameter *parameter1, BCParameter *parameter2)
TH1D *	MarginalizeByMetro (BCParameter *parameter)
virtual void	MCMCUserIterationInterface ()
void	ResetVarlist (int v)
double	SAHelperGetRadialCauchy ()
std::vector< double >	SAHelperGetRandomPointOnHypersphere ()
double	SAHelperSinusToNIntegral (int dim, double theta)
void	SAInitialize ()
double	SATemperature (double t)
double	SATemperatureBoltzmann (double t)
double	SATemperatureCauchy (double t)
virtual double	SATemperatureCustom (double t)
void	UnsetVar (int index)
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;
   fRandom = new TRandom3(0);

   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 84 of file BCIntegrate.cxx.

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

   fNbins=100;

   this->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 126 of file BCIntegrate.cxx.

{
   DeleteVarList();

   fx=0;

   delete fRandom;
   fRandom=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 1863 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::Out(BCLog::error,BCLog::error,"!!! This version of BAT is compiled without Cuba.");
   BCLog::Out(BCLog::error,BCLog::error,"    Use other integration methods or install Cuba and recompile BAT.");
#endif
}

double BCIntegrate::CubaIntegrate

(

)

Definition at line 1901 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::Out(BCLog::detail, BCLog::detail, Form(" --> Running Cuba/Vegas integation over %i dimensions.", fNvar));
   BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Maximum number of iterations: %i", (int)MAXEVAL));
   BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Aimed relative precision:     %e", EPSREL));

   return this -> CubaIntegrate(0, parameters_double, parameters_int);
#else
   BCLog::Out(BCLog::error,BCLog::error,"!!! This version of BAT is compiled without Cuba.");
   BCLog::Out(BCLog::error,BCLog::error,"    Use other integration methods or install Cuba and recompile BAT.");
   return 0.;
#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 1938 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::Out(BCLog::error,BCLog::error,"!!! This version of BAT is compiled without Cuba.");
   BCLog::Out(BCLog::error,BCLog::error,"    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 207 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 306 of file BCIntegrate.cxx.

{
   BCLog::Out(BCLog::warning, BCLog::warning, "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 333 of file BCIntegrate.cxx.

{
   double val=this->Eval(x);
   BCLog::Out(BCLog::detail, BCLog::detail, 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 320 of file BCIntegrate.cxx.

{
   BCLog::Out(BCLog::warning, BCLog::warning, "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 1802 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

Definition at line 1821 of file BCIntegrate.cxx.

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

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

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

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

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

         fBestFitParameters.clear();
         fBestFitParameterErrors.clear();
         fBestFitParameters = this -> MCMCGetMaximumPoint(i);

         for (int j = 0; j < fNvar; ++j)
            fBestFitParameterErrors.push_back(0.0); 

         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 1304 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.0,
                                 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 = this -> 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 = this -> 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 1410 of file BCIntegrate.cxx.

{
   // note: if f(x) is the function to be minimized, then
   // f(x) := - this->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 = this->LogEval(x);

   // run while still "hot enough"
   while ( this->SATemperature(t) > fSATmin )
   {
      // generate new state
      y = this->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 = this->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)
                  / this->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 = this -> 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 = this -> Eval( fBestFitParameters ); 
      }

   double probmaxsa = this -> 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(0.0); 
            }

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

   return;
}

virtual double BCIntegrate::FitFunction	(	std::vector< double >	x,
		std::vector< double >	parameters
	)			[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 490 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 1008 of file BCIntegrate.cxx.

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

   if(parIndex<0 || parIndex>fNvar)
   {
      BCLog::Out(BCLog::warning, BCLog::warning,
         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 1074 of file BCIntegrate.cxx.

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

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

   if(hindex>(int)fHProb2D.size()-1)
   {
      BCLog::Out(BCLog::warning, BCLog::warning,
         "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 1028 of file BCIntegrate.cxx.

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

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

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

   if(parIndex1>parIndex2)
   {
      BCLog::Out(BCLog::warning, BCLog::warning,
         "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::Out(BCLog::warning, BCLog::warning,
      "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 1294 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 1587 of file BCIntegrate.cxx.

{
   // do we have Cauchy (default), Boltzmann or custom annealing schedule?
   if (this->fSASchedule == BCIntegrate::kSABoltzmann)
      return this->GetProposalPointSABoltzmann(x, t);
   else if (this->fSASchedule == BCIntegrate::kSACauchy)
      return this->GetProposalPointSACauchy(x, t);
   else
      return this->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 1600 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]) * this->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 1618 of file BCIntegrate.cxx.

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

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

      norm = (fMax[0] - fMin[0]) * this->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 = this->SAHelperGetRandomPointOnHypersphere();

      // scale the vector by a random factor determined by the radial
      // part of the fNvar-dimensional Cauchy distribution
      double radial = this->SATemperature(t) * this->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 1655 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 1098 of file BCIntegrate.cxx.

{
   this->GetRandomVector(x);

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

   return this->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 1109 of file BCIntegrate.cxx.

{
   double p = 1.1;
   double g = 1.0;

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

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

      p = fRandom->Rndm();

      g = this -> EvalSampling(x);
   }

   return this->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 1154 of file BCIntegrate.cxx.

{
   // get new point
   this->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 = this -> 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 1215 of file BCIntegrate.cxx.

{
   // get new point
   this->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 = this -> LogEvalSampling(fXmetro0); // old point
   double p1=0; // new point
   if(in)
      p1= this -> 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 1266 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 1280 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; };

double* BCIntegrate::GetSAP2LogProb

(

)

[inline]

Definition at line 229 of file BCIntegrate.h.

         { return &fSALogProb; }; 

int* BCIntegrate::GetSAP2NIterations

(

)

[inline]

Definition at line 235 of file BCIntegrate.h.

         { return &fSANIterations; };

double* BCIntegrate::GetSAP2Temperature

(

)

[inline]

Definition at line 226 of file BCIntegrate.h.

         { return &fSATemperature; };

std::vector<double>* BCIntegrate::GetSAP2x

(

)

[inline]

Definition at line 232 of file BCIntegrate.h.

         { return &fSAx; }; 

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; }

void BCIntegrate::InitMetro

(

)

Initializes the Metropolis algorithm (for details see manual)

Definition at line 1130 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 = this ->  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 560 of file BCIntegrate.cxx.

{
   // print debug information
   BCLog::Out(BCLog::debug, BCLog::debug, 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::Out(BCLog::detail, BCLog::detail, 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
      this -> GetRandomPointImportance(randx);

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

      // calculate sampling distributions at that point
      double val_g = this -> 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::Out(BCLog::debug, BCLog::debug,
            Form("BCIntegrate::IntegralImportance. Iteration %d, integral: %d.", i+1, sumI/(i+1)));
   }

   // calculate integral
   double result=sumI/Niter;

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

   return result;
}

double BCIntegrate::IntegralMC

(

std::vector< double >

x

)

Definition at line 393 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::Out(BCLog::detail, BCLog::detail, Form(" --> Running MC integation over %i dimensions.", NvarNow));
   BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Maximum number of iterations: %i", this->GetNIterationsMax()));
   BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Aimed relative precision:     %e", this->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
      this -> 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 = this->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::Out(BCLog::detail, BCLog::detail, Form(" --> Result of integration:        %e +- %e   in %i iterations.", integral, sqrt(variance), fNIterations));
   BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Obtained relative precision:  %e. ", sqrt(variance) / integral));

   return integral;
}

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 386 of file BCIntegrate.cxx.

{
   this->SetVarList(varlist);
   return IntegralMC(x);
}

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 491 of file BCIntegrate.cxx.

{
   // print debug information
   BCLog::Out(BCLog::debug, BCLog::debug, 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::Out(BCLog::detail, BCLog::detail, 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
      this -> GetRandomPointSamplingMetro(randx);

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

      // calculate sampling distributions at that point
      double val_g = this -> 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::Out(BCLog::debug, BCLog::debug,
            Form("BCIntegrate::IntegralMetro. Iteration %d, integral: %d.", i+1, sumI/(i+1)));
   }

   // calculate integral
   double result=sumI/Niter;

   // print debug information
   BCLog::Out(BCLog::detail, BCLog::detail, 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 354 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 this -> IntegralMetro(parameter);

      case BCIntegrate::kIntImportance:
         return this -> IntegralImportance(parameter);

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

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

   return 0;
}

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 313 of file BCIntegrate.cxx.

{
   // this method should better also be overloaded
   return log(this->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 327 of file BCIntegrate.cxx.

{
   return log(this->EvalSampling(x));
}

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 648 of file BCIntegrate.cxx.

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

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

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

   return 0;
}

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 627 of file BCIntegrate.cxx.

{
   BCLog::Out(BCLog::detail, BCLog::detail,
                   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::Out(BCLog::warning, BCLog::warning,
      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 816 of file BCIntegrate.cxx.

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

   BCLog::Out(BCLog::detail, BCLog::detail,
         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::Out(BCLog::detail, BCLog::detail,
         Form(" --> Marginalizing %d 1D distributions and %d 2D distributions.", fNvar, nh2d));

   // prepare function fitting
   double dx = 0.0;
   double dy = 0.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::Out(BCLog::debug, BCLog::debug,
            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 = this -> 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 = this -> 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::Out(BCLog::detail, BCLog::detail,
         Form("BCIntegrate::MarginalizeAllByMetro done with %i trials and %i accepted trials. Efficiency is %f",fNmetro, fNacceptedMCMC, double(fNacceptedMCMC)/double(fNmetro)));

   return fNvar+nh2d;
}

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 698 of file BCIntegrate.cxx.

{
   // set parameter to marginalize
   this->ResetVarlist(1);
   int index1 = parameter1->GetIndex();
   this->UnsetVar(index1);
   int index2 = parameter2->GetIndex();
   this->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::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 666 of file BCIntegrate.cxx.

{
   // set parameter to marginalize
   this->ResetVarlist(1);
   int index = parameter->GetIndex();
   this->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.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::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 779 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;
}

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 743 of file BCIntegrate.cxx.

{
   int niter = fMarkovChainNIterations;

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

   BCLog::Out(BCLog::detail, BCLog::detail,
      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;
}

void BCIntegrate::MCMCFillErrorBand

(

)

[private]

Definition at line 2056 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 < this -> MCMCGetNChains(); ++ichain)
         {
            // calculate y
            double y = this -> FitFunction(xvec, this -> 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 < this -> MCMCGetNChains(); ++ichain)
         {
            // calculate y
            double y = this -> FitFunction(xvec, this -> MCMCGetx(ichain));

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

         xvec.clear();
      }
   }
}

void BCIntegrate::MCMCIterationInterface

(

)

[private, virtual]

Reimplemented from BCEngineMCMC.

Definition at line 2043 of file BCIntegrate.cxx.

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

   // fill error band
   this -> MCMCFillErrorBand();

   // do user defined stuff
   this -> MCMCUserIterationInterface();
}

virtual void BCIntegrate::MCMCUserIterationInterface

(

)

[inline, virtual]

Reimplemented in BCGoFTest.

Definition at line 718 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 299 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 1711 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 = this->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] = this->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 1663 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 1772 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
         * this->SAHelperSinusToNIntegral(dim - 2, theta);
}

void BCIntegrate::SAInitialize

(

)

Definition at line 2122 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 1552 of file BCIntegrate.cxx.

{
   // do we have Cauchy (default), Boltzmann or custom annealing schedule?
   if (this->fSASchedule == BCIntegrate::kSABoltzmann)
      return this->SATemperatureBoltzmann(t);
   else if (this->fSASchedule == BCIntegrate::kSACauchy)
      return this->SATemperatureCauchy(t);
   else
      return this->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 1565 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 1572 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 1579 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 350 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 362 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 356 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 368 of file BCIntegrate.h.

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

void BCIntegrate::SetErrorBandHisto

(

TH2D *

h

)

[inline]

Sets errorband histogram

Definition at line 408 of file BCIntegrate.h.

         { fErrorBandXY = h; };

void BCIntegrate::SetFillErrorBand

(

bool

flag = true

)

[inline]

Definition at line 322 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 334 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 340 of file BCIntegrate.h.

         { fFitFunctionIndexY = index; };

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

Definition at line 343 of file BCIntegrate.h.

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

void BCIntegrate::SetFlagIgnorePrevOptimization

(

bool

flag

)

[inline]

Definition at line 247 of file BCIntegrate.h.

      { fFlagIgnorePrevOptimization = flag; }; 

void BCIntegrate::SetFlagWriteSAToFile

(

bool

flag

)

[inline]

Definition at line 421 of file BCIntegrate.h.

      { fFlagWriteSAToFile = flag; }; 

void BCIntegrate::SetIntegrationMethod

(

BCIntegrate::BCIntegrationMethod

method

)

Parameters:

method

The integration method

Definition at line 342 of file BCIntegrate.cxx.

{
#ifdef HAVE_CUBA_H
   fIntegrationMethod = method;
#else
   BCLog::Out(BCLog::warning,BCLog::warning,"!!! This version of BAT is compiled without Cuba.");
   BCLog::Out(BCLog::warning,BCLog::warning,"    Monte Carlo Sampled Mean integration method will be used.");
   BCLog::Out(BCLog::warning,BCLog::warning,"    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 269 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 399 of file BCIntegrate.h.

         { fMarkovChainAutoN = flag; };

void BCIntegrate::SetMarkovChainInitialPosition

(

std::vector< double >

position

)

Sets the initial position for the Markov chain

Definition at line 196 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 393 of file BCIntegrate.h.

         { fMarkovChainNIterations = niterations;
           fMarkovChainAutoN = false; }

void BCIntegrate::SetMarkovChainStepSize

(

double

stepsize

)

[inline]

Sets the step size for Markov chains

Definition at line 388 of file BCIntegrate.h.

         { fMarkovChainStepSize = stepsize; };

void BCIntegrate::SetMarkovChainTree

(

TTree *

tree

)

[inline]

Sets the ROOT tree containing the Markov chain

Definition at line 379 of file BCIntegrate.h.

         { fMarkovChainTree = tree; };

void BCIntegrate::SetMinuitArlist

(

double *

arglist

)

[inline]

Definition at line 243 of file BCIntegrate.h.

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

void BCIntegrate::SetMode

(

std::vector< double >

mode

)

Sets mode

Definition at line 1790 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 226 of file BCIntegrate.cxx.

{
   if (fNvar == 0)
      return;

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

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

   fMCMCH1NBins[index] = nbins;

   return;

//    if(n<2)
//    {
//       BCLog::Out(BCLog::warning, BCLog::warning,"BCIntegrate::SetNbins. Number of bins less than 2 makes no sense. Setting to 2.");
//       n=2;
//    }
//    this -> 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 295 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 284 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 290 of file BCIntegrate.h.

         { fNSamplesPer2DBin = n; };

void BCIntegrate::SetOptimizationMethod

(

BCIntegrate::BCOptimizationMethod

method

)

[inline]

Parameters:

method

The mode finding method

Definition at line 274 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 154 of file BCIntegrate.cxx.

{
   DeleteVarList();

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

   this->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 301 of file BCIntegrate.h.

         { fRelativePrecision = relprecision; };

void BCIntegrate::SetSASchedule

(

BCIntegrate::BCSASchedule

schedule

)

[inline]

Parameters:

method

The Simulated Annealing schedule

Definition at line 279 of file BCIntegrate.h.

         { fSASchedule = schedule; };

void BCIntegrate::SetSAT0

(

double

T0

)

[inline]

Parameters:

T0

new value for Simulated Annealing starting temperature.

Definition at line 413 of file BCIntegrate.h.

         { fSAT0 = T0; }

void BCIntegrate::SetSATmin

(

double

Tmin

)

[inline]

Parameters:

Tmin

new value for Simulated Annealing threshold temperature.

Definition at line 418 of file BCIntegrate.h.

         { fSATmin = Tmin; }

void BCIntegrate::SetSATree

(

TTree *

tree

)

[inline]

Definition at line 426 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 261 of file BCIntegrate.h.

{fVarlist[index]=1;};

void BCIntegrate::SetVarList

(

int *

varlist

)

Parameters:

varlist

A list of parameters

Definition at line 292 of file BCIntegrate.cxx.

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

void BCIntegrate::UnsetFillErrorBand

(

)

[inline]

Definition at line 328 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 446 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 374 of file BCIntegrate.h.

         { fFlagWriteMarkovChain = flag; fMCMCFlagWriteChainToFile = 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 838 of file BCIntegrate.h.

BCDataPoint* BCIntegrate::fDataPointLowerBoundaries [protected]

data point containing the lower boundaries of possible data values

Definition at line 832 of file BCIntegrate.h.

BCDataPoint* BCIntegrate::fDataPointUpperBoundaries [protected]

data point containing the upper boundaries of possible data values

Definition at line 836 of file BCIntegrate.h.

double BCIntegrate::fError [private]

The uncertainty in the most recent Monte Carlo integration

Definition at line 780 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 747 of file BCIntegrate.h.

BCIntegrate::BCMarginalizationMethod BCIntegrate::fMarginalizationMethod [private]

Current marginalization method

Definition at line 751 of file BCIntegrate.h.

bool BCIntegrate::fMarkovChainAutoN [protected]

Definition at line 828 of file BCIntegrate.h.

int BCIntegrate::fMarkovChainNIterations [protected]

Definition at line 826 of file BCIntegrate.h.

double BCIntegrate::fMarkovChainStepSize [protected]

Step size in the Markov chain relative to min and max

Definition at line 824 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 797 of file BCIntegrate.h.

double* BCIntegrate::fMax [private]

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

Definition at line 735 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 731 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 785 of file BCIntegrate.h.

int BCIntegrate::fNbins [protected]

Number of bins per dimension for the marginalized distributions

Definition at line 815 of file BCIntegrate.h.

int BCIntegrate::fNIterations [private]

Number of iterations in the most recent Monte Carlo integation

Definition at line 772 of file BCIntegrate.h.

int BCIntegrate::fNIterationsMax [private]

Maximum number of iterations

Definition at line 768 of file BCIntegrate.h.

int BCIntegrate::fNiterPerDimension [private]

Number of iteration per dimension for Monte Carlo integration.

Definition at line 743 of file BCIntegrate.h.

int BCIntegrate::fNmetro [private]

The number of iterations in the Metropolis integration

Definition at line 784 of file BCIntegrate.h.

int BCIntegrate::fNSamplesPer2DBin [protected]

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

Definition at line 820 of file BCIntegrate.h.

int BCIntegrate::fNvar [protected]

Number of variables to integrate over.

Definition at line 811 of file BCIntegrate.h.

BCIntegrate::BCOptimizationMethod BCIntegrate::fOptimizationMethod [private]

Current mode finding method

Definition at line 755 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 760 of file BCIntegrate.h.

TRandom3* BCIntegrate::fRandom [protected]

A ROOT random number generator

Definition at line 842 of file BCIntegrate.h.

double BCIntegrate::fRelativePrecision [private]

Relative precision aimed at in the Monte Carlo integation

Definition at line 776 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 764 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 739 of file BCIntegrate.h.

BCParameterSet* BCIntegrate::fx [private]

Set of parameters for the integration.

Definition at line 719 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 789 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 793 of file BCIntegrate.h.

---

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

• BCIntegrate.h
• BCIntegrate.cxx