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

#include <BCIntegrate.h>

Inheritance diagram for BCIntegrate:

Collaboration diagram for BCIntegrate:

Public Types
Enumerators
enum	BCOptimizationMethod { kOptEmpty, kOptSimAnn, kOptMetropolis, kOptMinuit, kOptDefault, NOptMethods }

enum	BCIntegrationMethod { kIntEmpty, kIntMonteCarlo, kIntCuba, kIntGrid, kIntDefault, NIntMethods }

enum	BCMarginalizationMethod { kMargEmpty, kMargMetropolis, kMargMonteCarlo, kMargGrid, kMargDefault, NMargMethods }

enum	BCSASchedule { kSACauchy, kSABoltzmann, kSACustom, NSAMethods }

enum	BCCubaMethod { kCubaVegas, kCubaSuave, kCubaDivonne, kCubaCuhre, NCubaMethods }

Function pointer types
typedef void(BCIntegrate::*	tRandomizer )(std::vector< double > &) const

typedef double(BCIntegrate::*	tEvaluator )(std::vector< double > &, const std::vector< double > &, bool &)

typedef void(*	tIntegralUpdater )(const std::vector< double > &, const int &, double &, double &)

Public Types inherited from BCEngineMCMC
enum	Precision { kLow, kMedium, kHigh, kVeryHigh }

Public Member Functions
Constructors and destructors
	BCIntegrate ()

	BCIntegrate (const BCIntegrate &bcintegrate)

virtual	~BCIntegrate ()

Assignment operators
BCIntegrate &	operator= (const BCIntegrate &bcintegrate)

Member functions (get)
int	ReadMarginalizedFromFile (const char *file)

BCH1D *	GetMarginalized (const BCParameter *parameter)

BCH1D *	GetMarginalized (const char *name)

BCH1D *	GetMarginalized (unsigned index)

BCH2D *	GetMarginalized (const BCParameter parameter1, const BCParameter parameter2)

BCH2D *	GetMarginalized (const char name1, const char name2)

BCH2D *	GetMarginalized (unsigned index1, unsigned index2)

int	PrintAllMarginalized1D (const char *filebase)

int	PrintAllMarginalized2D (const char *filebase)

int	PrintAllMarginalized (const char *file, std::string options1d="BTsiB3CS1D0pdf0Lmeanmode", std::string options2d="BTfB3CS1meangmode", unsigned int hdiv=1, unsigned int ndiv=1)

double	GetIntegral () const

BCIntegrate::BCOptimizationMethod	GetOptimizationMethod () const

BCIntegrate::BCIntegrationMethod	GetIntegrationMethod () const

BCIntegrate::BCMarginalizationMethod	GetMarginalizationMethod () const

BCIntegrate::BCSASchedule	GetSASchedule () const

void	GetRandomVectorUnitHypercube (std::vector< double > &x) const

void	GetRandomVectorInParameterSpace (std::vector< double > &x) const

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

int	GetNIterationsMin () const

int	GetNIterationsMax () const

int	GetNIterationsPrecisionCheck () const

int	GetNIterationsOutput () const

int	GetNIterations () const

double	GetRelativePrecision () const

double	GetAbsolutePrecision () const

BCCubaMethod	GetCubaIntegrationMethod () const

const BCCubaOptions::Vegas &	GetCubaVegasOptions () const

const BCCubaOptions::Suave &	GetCubaSuaveOptions () const

const BCCubaOptions::Divonne &	GetCubaDivonneOptions () const

const BCCubaOptions::Cuhre &	GetCubaCuhreOptions () const

BCH1D *	GetSlice (const BCParameter *parameter, const std::vector< double > parameters=std::vector< double >(0), int bins=0, bool flag_norm=true)

BCH1D *	GetSlice (const char *name, const std::vector< double > parameters=std::vector< double >(0), int nbins=0, bool flag_norm=true)

BCH2D *	GetSlice (const BCParameter parameter1, const BCParameter parameter2, const std::vector< double > parameters=std::vector< double >(0), int bins=0, bool flag_norm=true)

BCH2D *	GetSlice (const char name1, const char name2, const std::vector< double > parameters=std::vector< double >(0), int nbins=0, bool flag_norm=true)

BCH2D *	GetSlice (unsigned index1, unsigned index2, const std::vector< double > parameters=std::vector< double >(0), int nbins=0, bool flag_norm=true)

double	GetError () const

TMinuit *	GetMinuit ()

int	GetMinuitErrorFlag () const

double	GetSAT0 () const

double	GetSATmin () const

double	GetBestFitParameter (unsigned index) const

double	GetBestFitParameterError (unsigned index) const

double	GetLogMaximum ()

const std::vector< double > &	GetBestFitParameters () const

const std::vector< double > &	GetBestFitParameterErrors () const

Member functions (set)
void	SetMinuitArlist (double *arglist)

void	SetFlagIgnorePrevOptimization (bool flag)

void	SetOptimizationMethod (BCIntegrate::BCOptimizationMethod method)

void	SetIntegrationMethod (BCIntegrate::BCIntegrationMethod method)

void	SetMarginalizationMethod (BCIntegrate::BCMarginalizationMethod method)

void	SetSASchedule (BCIntegrate::BCSASchedule schedule)

void	SetNIterationsMin (int niterations)

void	SetNIterationsMax (int niterations)

void	SetNIterationsPrecisionCheck (int niterations)

void	SetNIterationsOutput (int niterations)

void	SetRelativePrecision (double relprecision)

void	SetAbsolutePrecision (double absprecision)

void	SetCubaIntegrationMethod (BCCubaMethod type)

void	SetCubaOptions (const BCCubaOptions::Vegas &options)

void	SetCubaOptions (const BCCubaOptions::Suave &options)

void	SetCubaOptions (const BCCubaOptions::Divonne &options)

void	SetCubaOptions (const BCCubaOptions::Cuhre &options)

void	SetSAT0 (double T0)

void	SetSATmin (double Tmin)

void	SetFlagWriteSAToFile (bool flag)

TTree *	GetSATree ()

void	InitializeSATree ()

Public Member Functions inherited from BCEngineMCMC
void	WriteMarkovChain (bool flag)

	BCEngineMCMC ()

	BCEngineMCMC (const BCEngineMCMC &enginemcmc)

virtual	~BCEngineMCMC ()

BCEngineMCMC &	operator= (const BCEngineMCMC &engineMCMC)

unsigned	MCMCGetNChains () const

unsigned	MCMCGetNLag () const

const std::vector< unsigned > &	MCMCGetNIterations () const

unsigned	MCMCGetCurrentIteration () const

unsigned	MCMCGetCurrentChain () const

unsigned	MCMCGetNIterationsConvergenceGlobal () const

bool	MCMCGetFlagConvergenceGlobal () const

unsigned	MCMCGetNIterationsMax () const

unsigned	MCMCGetNIterationsRun () const

unsigned	MCMCGetNIterationsPreRunMin () const

unsigned	MCMCGetNIterationsUpdate () const

unsigned	MCMCGetNIterationsUpdateMax () const

std::vector< int >	MCMCGetNTrialsTrue () const

int	MCMCGetNTrials () const

const std::vector< double > &	MCMCGetprobMean () const

const std::vector< double > &	MCMCGetVariance () const

const std::vector< double > &	MCMCGetTrialFunctionScaleFactor () const

std::vector< double >	MCMCGetTrialFunctionScaleFactor (unsigned ichain) const

double	MCMCGetTrialFunctionScaleFactor (unsigned ichain, unsigned ipar)

const std::vector< double > &	MCMCGetx () const

std::vector< double >	MCMCGetx (unsigned ichain)

double	MCMCGetx (unsigned ichain, unsigned ipar) const

const std::vector< double > &	MCMCGetLogProbx () const

double	MCMCGetLogProbx (unsigned ichain)

int	MCMCGetPhase () const

const std::vector< double > &	MCMCGetMaximumPoints () const

std::vector< double >	MCMCGetMaximumPoint (unsigned i) const

const std::vector< double > &	MCMCGetMaximumLogProb () const

int	MCMCGetFlagInitialPosition () const

double	MCMCGetRValueCriterion () const

double	MCMCGetRValueParametersCriterion () const

double	MCMCGetRValue () const

double	MCMCGetRValueParameters (unsigned i)

bool	MCMCGetRValueStrict () const

bool	MCMCGetFlagRun () const

TTree *	MCMCGetMarkovChainTree (unsigned i)

BCH1D *	MCMCGetH1Marginalized (unsigned i)

BCH2D *	MCMCGetH2Marginalized (unsigned i, unsigned j)

BCParameter *	GetParameter (int index) const

BCParameter *	GetParameter (const char *name) const

unsigned int	GetNParameters () const

unsigned int	GetNFixedParameters ()

unsigned int	GetNFreeParameters ()

const std::vector< double > &	GetBestFitParametersMarginalized () const

void	MCMCSetTrialFunctionScaleFactor (std::vector< double > scale)

void	MCMCSetNChains (unsigned n)

void	MCMCSetNLag (unsigned n)

void	MCMCSetNIterationsMax (unsigned n)

void	MCMCSetNIterationsRun (unsigned n)

void	MCMCSetNIterationsPreRunMin (unsigned n)

void	MCMCSetNIterationsUpdate (unsigned n)

void	MCMCSetNIterationsUpdateMax (unsigned n)

void	MCMCSetMinimumEfficiency (double efficiency)

void	MCMCSetMaximumEfficiency (double efficiency)

void	MCMCSetRandomSeed (unsigned seed)

void	MCMCSetWriteChainToFile (bool flag)

void	MCMCSetWritePreRunToFile (bool flag)

void	MCMCSetInitialPositions (const std::vector< double > &x0s)

void	MCMCSetInitialPositions (std::vector< std::vector< double > > x0s)

void	MCMCSetFlagInitialPosition (int flag)

void	MCMCSetFlagOrderParameters (bool flag)

void	MCMCSetFlagFillHistograms (bool flag)

void	MCMCSetFlagPreRun (bool flag)

void	MCMCSetRValueCriterion (double r)

void	MCMCSetRValueParametersCriterion (double r)

void	MCMCSetRValueStrict (bool strict=true)

void	MCMCSetMarkovChainTrees (const std::vector< TTree * > &trees)

void	MCMCInitializeMarkovChainTrees ()

int	SetMarginalized (unsigned index, TH1D *h)

int	SetMarginalized (unsigned index1, unsigned index2, TH2D *h)

void	MCMCSetValuesDefault ()

void	MCMCSetValuesQuick ()

void	MCMCSetValuesDetail ()

void	MCMCSetPrecision (BCEngineMCMC::Precision precision)

void	SetNbins (unsigned int nbins)

void	Copy (const BCEngineMCMC &enginemcmc)

virtual int	AddParameter (const char name, double min, double max, const char latexname="")

virtual int	AddParameter (BCParameter *parameter)

virtual void	MCMCTrialFunction (unsigned ichain, std::vector< double > &x)

virtual double	MCMCTrialFunctionSingle (unsigned ichain, unsigned ipar)

bool	MCMCGetProposalPointMetropolis (unsigned chain, std::vector< double > &x)

bool	MCMCGetProposalPointMetropolis (unsigned chain, unsigned parameter, std::vector< double > &x)

bool	MCMCGetNewPointMetropolis (unsigned chain=0)

bool	MCMCGetNewPointMetropolis (unsigned chain, unsigned parameter)

void	MCMCInChainCheckMaximum ()

void	MCMCInChainUpdateStatistics ()

void	MCMCInChainFillHistograms ()

void	MCMCInChainTestConvergenceAllChains ()

void	MCMCInChainWriteChains ()

int	MCMCMetropolis ()

int	MCMCMetropolisPreRun ()

void	MCMCResetRunStatistics ()

void	MCMCInitializeMarkovChains ()

int	MCMCInitialize ()

void	ClearParameters (bool hard=false)

virtual void	MCMCIterationInterface ()

virtual void	MCMCUserIterationInterface ()

virtual void	MCMCCurrentPointInterface (std::vector< double > &, int, bool)

Protected Member Functions
unsigned	IntegrationOutputFrequency () const

void	LogOutputAtStartOfIntegration (BCIntegrationMethod type, BCCubaMethod cubatype)

void	LogOutputAtIntegrationStatusUpdate (BCIntegrationMethod type, double integral, double absprecision, int nIterations)

void	LogOutputAtEndOfIntegration (double integral, double absprecision, double relprecision, int nIterations)

void	Copy (const BCIntegrate &bcintegrate)

Protected Attributes
int	fID

TMinuit *	fMinuit

double	fMinuitArglist [2]

int	fMinuitErrorFlag

bool	fFlagIgnorePrevOptimization

double	fSAT0

double	fSATmin

TTree *	fTreeSA

bool	fFlagWriteSAToFile

int	fSANIterations

double	fSATemperature

double	fSALogProb

std::vector< double >	fSAx

std::vector< BCH1D * >	fMarginalized1D

std::vector< BCH2D * >	fMarginalized2D

bool	fFlagMarginalized

Protected Attributes inherited from BCEngineMCMC
BCParameterSet	fParameters

unsigned	fMCMCNChains

unsigned	fMCMCNLag

std::vector< unsigned >	fMCMCNIterations

int	fMCMCCurrentIteration

int	fMCMCCurrentChain

unsigned	fMCMCNIterationsUpdate

unsigned	fMCMCNIterationsUpdateMax

int	fMCMCNIterationsConvergenceGlobal

bool	fMCMCFlagConvergenceGlobal

unsigned	fMCMCNIterationsMax

unsigned	fMCMCNIterationsRun

unsigned	fMCMCNIterationsPreRunMin

std::vector< int >	fMCMCNTrialsTrue

unsigned	fMCMCNTrials

bool	fMCMCFlagWriteChainToFile

bool	fMCMCFlagWritePreRunToFile

std::vector< double >	fMCMCTrialFunctionScaleFactor

std::vector< double >	fMCMCTrialFunctionScaleFactorStart

bool	fMCMCFlagPreRun

bool	fMCMCFlagRun

std::vector< double >	fMCMCInitialPosition

double	fMCMCEfficiencyMin

double	fMCMCEfficiencyMax

int	fMCMCFlagInitialPosition

bool	fMCMCFlagOrderParameters

int	fMCMCPhase

std::vector< double >	fMCMCx

std::vector< double >	fMCMCxMax

std::vector< double >	fMCMCxMean

std::vector< double >	fMCMCxVar

std::vector< double >	fMCMCprob

std::vector< double >	fMCMCprobMax

std::vector< double >	fMCMCprobMean

std::vector< double >	fMCMCprobVar

bool	fMCMCRValueUseStrict

double	fMCMCRValueCriterion

double	fMCMCRValueParametersCriterion

double	fMCMCRValue

std::vector< double >	fMCMCRValueParameters

TRandom3 *	fRandom

std::vector< TH1D * >	fMCMCH1Marginalized

std::vector< TH2D * >	fMCMCH2Marginalized

std::vector< TTree * >	fMCMCTrees

std::vector< double >	fMCMCBestFitParameters

double	fMCMCLogMaximum

std::vector< double >	fMarginalModes

Private Member Functions
std::vector< double >	FindModeMinuit (std::vector< double > &mode, std::vector< double > &errors, std::vector< double > start=std::vector< double >(0), int printlevel=1)

std::vector< double >	FindModeMCMC (std::vector< double > &mode, std::vector< double > &errors)

std::vector< double >	FindModeSA (std::vector< double > &mode, std::vector< double > &errors, std::vector< double > start=std::vector< double >(0))

double	IntegrateCuba ()

double	IntegrateCuba (BCCubaMethod cubatype)

double	IntegrateSlice ()

Private Attributes
BCIntegrate::BCOptimizationMethod	fOptimizationMethodCurrent

BCIntegrate::BCOptimizationMethod	fOptimizationMethodUsed

BCIntegrate::BCIntegrationMethod	fIntegrationMethodCurrent

BCIntegrate::BCIntegrationMethod	fIntegrationMethodUsed

BCIntegrate::BCMarginalizationMethod	fMarginalizationMethodCurrent

BCIntegrate::BCMarginalizationMethod	fMarginalizationMethodUsed

BCIntegrate::BCSASchedule	fSASchedule

unsigned	fNIterationsMin

unsigned	fNIterationsMax

unsigned	fNIterationsPrecisionCheck

unsigned	fNIterationsOutput

int	fNIterations

std::vector< double >	fBestFitParameters

std::vector< double >	fBestFitParameterErrors

double	fLogMaximum

double	fIntegral

double	fRelativePrecision

double	fAbsolutePrecision

double	fError

BCCubaMethod	fCubaIntegrationMethod

BCCubaOptions::Vegas	fCubaVegasOptions

BCCubaOptions::Suave	fCubaSuaveOptions

BCCubaOptions::Divonne	fCubaDivonneOptions

BCCubaOptions::Cuhre	fCubaCuhreOptions

Member functions (miscellaneous methods)
virtual double	Eval (const std::vector< double > &x)

virtual double	LogEval (const std::vector< double > &x)

double	Normalize ()

double	Integrate (BCIntegrationMethod intmethod)

double	Integrate ()

double	Integrate (BCIntegrationMethod type, tRandomizer randomizer, tEvaluator evaluator, tIntegralUpdater updater, std::vector< double > &sums)

double	EvaluatorMC (std::vector< double > &sums, const std::vector< double > &point, bool &accepted)

int	MarginalizeAll ()

int	MarginalizeAll (BCMarginalizationMethod margmethod)

virtual void	MarginalizePreprocess ()

virtual void	MarginalizePostprocess ()

void	SAInitialize ()

std::vector< double >	FindMode (std::vector< double > start=std::vector< double >())

std::vector< double >	FindMode (BCIntegrate::BCOptimizationMethod optmethod, std::vector< double > start=std::vector< double >())

double	SATemperature (double t)

double	SATemperatureBoltzmann (double t)

double	SATemperatureCauchy (double t)

virtual double	SATemperatureCustom (double t)

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

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

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

virtual std::vector< double >	GetProposalPointSACustom (const std::vector< double > &x, int t)

std::vector< double >	SAHelperGetRandomPointOnHypersphere ()

double	SAHelperGetRadialCauchy ()

double	SAHelperSinusToNIntegral (int dim, double theta)

virtual void	ResetResults ()

std::string	DumpIntegrationMethod (BCIntegrationMethod type)

std::string	DumpCurrentIntegrationMethod ()

std::string	DumpUsedIntegrationMethod ()

std::string	DumpMarginalizationMethod (BCMarginalizationMethod type)

std::string	DumpCurrentMarginalizationMethod ()

std::string	DumpUsedMarginalizationMethod ()

std::string	DumpOptimizationMethod (BCOptimizationMethod type)

std::string	DumpCurrentOptimizationMethod ()

std::string	DumpUsedOptimizationMethod ()

std::string	DumpCubaIntegrationMethod (BCCubaMethod type)

std::string	DumpCubaIntegrationMethod ()

void	SetBestFitParameters (const std::vector< double > &x)

void	SetBestFitParameters (const std::vector< double > &x, const double &new_value, double &old_value)

unsigned	GetNIntegrationVariables ()

double	CalculateIntegrationVolume ()

bool	CheckMarginalizationAvailability (BCMarginalizationMethod type)

bool	CheckMarginalizationIndices (TH1 *hist, const std::vector< unsigned > &index)

static void	IntegralUpdaterMC (const std::vector< double > &sums, const int &nIterations, double &integral, double &absprecision)

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

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

Member Typedef Documentation

typedef double(BCIntegrate::* BCIntegrate::tEvaluator)(std::vector< double > &, const std::vector< double > &, bool &)

A pointer for a function that evaluates at a point

Definition at line 141 of file BCIntegrate.h.

typedef void(* BCIntegrate::tIntegralUpdater)(const std::vector< double > &, const int &, double &, double &)

A pointer for a function that updates the integral and absolute precision

Definition at line 145 of file BCIntegrate.h.

typedef void(BCIntegrate::* BCIntegrate::tRandomizer)(std::vector< double > &) const

A pointer for a function that chooses a next random point

Definition at line 137 of file BCIntegrate.h.

Member Enumeration Documentation

enum BCIntegrate::BCCubaMethod

An enumerator for Cuba integration methods

Enumerator
kCubaVegas
kCubaSuave
kCubaDivonne
kCubaCuhre
NCubaMethods

Definition at line 128 of file BCIntegrate.h.

128 { kCubaVegas, kCubaSuave, kCubaDivonne, kCubaCuhre, NCubaMethods};

enum BCIntegrate::BCIntegrationMethod

An enumerator for integration algorithms

Enumerator
kIntEmpty
kIntMonteCarlo
kIntCuba
kIntGrid
kIntDefault
NIntMethods

Definition at line 104 of file BCIntegrate.h.

                             {
       kIntEmpty,
       kIntMonteCarlo,
       kIntCuba,
       kIntGrid,
       kIntDefault,
       NIntMethods };

enum BCIntegrate::BCMarginalizationMethod

An enumerator for marginalization algorithms

Enumerator
kMargEmpty
kMargMetropolis
kMargMonteCarlo
kMargGrid
kMargDefault
NMargMethods

Definition at line 114 of file BCIntegrate.h.

                                 {
      kMargEmpty,
      kMargMetropolis,
      kMargMonteCarlo,
      kMargGrid,
      kMargDefault,
      NMargMethods };

enum BCIntegrate::BCOptimizationMethod

An enumerator for the mode finding algorithm

Enumerator
kOptEmpty
kOptSimAnn
kOptMetropolis
kOptMinuit
kOptDefault
NOptMethods

Definition at line 94 of file BCIntegrate.h.

                             {
     kOptEmpty,
     kOptSimAnn,
     kOptMetropolis,
     kOptMinuit,
     kOptDefault,
     NOptMethods };

enum BCIntegrate::BCSASchedule

An enumerator for the Simulated Annealing schedule

Enumerator
kSACauchy
kSABoltzmann
kSACustom
NSAMethods

Definition at line 124 of file BCIntegrate.h.

124 { kSACauchy, kSABoltzmann, kSACustom, NSAMethods };

Constructor & Destructor Documentation

BCIntegrate::BCIntegrate ( )

A constructor

Definition at line 70 of file BCIntegrate.cxx.

                          : BCEngineMCMC(),
                              fID(BCLog::GetHIndex()),
                              fMinuit(0),
                              fMinuitErrorFlag(0),
                              fFlagIgnorePrevOptimization(false),
                              fSAT0(100),
                              fSATmin(0.1),
                              fTreeSA(0),
                              fFlagWriteSAToFile(false),
                              fSANIterations(0),
                              fSATemperature(0),
                              fSALogProb(0),
                              fMarginalized1D(0),
                              fMarginalized2D(0),
                              fFlagMarginalized(false),
                              fOptimizationMethodCurrent(BCIntegrate::kOptDefault),
                              fOptimizationMethodUsed(BCIntegrate::kOptEmpty),
                              fIntegrationMethodCurrent(BCIntegrate::kIntDefault),
                              fIntegrationMethodUsed(BCIntegrate::kIntEmpty),
                              fMarginalizationMethodCurrent(BCIntegrate::kMargDefault),
                              fMarginalizationMethodUsed(BCIntegrate::kMargEmpty),
                              fSASchedule(BCIntegrate::kSACauchy),
                              fNIterationsMin(0),
                              fNIterationsMax(1000000),
                              fNIterationsPrecisionCheck(1000),
                              fNIterationsOutput(0),
                              fNIterations(0),
                              fLogMaximum(-std::numeric_limits<double>::max()),
                              fIntegral(-1),
                              fRelativePrecision(1e-2),
                              fAbsolutePrecision(1e-6),
                              fError(-999.),
                              fCubaIntegrationMethod(BCIntegrate::kCubaVegas)
 {
    fMinuitArglist[0] = 20000;
    fMinuitArglist[1] = 0.01;
 }

BCIntegrate::BCIntegrate ( const BCIntegrate & bcintegrate)

The copy constructor

Definition at line 109 of file BCIntegrate.cxx.

                                                   : BCEngineMCMC(other)
 {
    Copy(other);
 }

BCIntegrate::~BCIntegrate ( )

virtual

The default destructor

Definition at line 170 of file BCIntegrate.cxx.

 {
    delete fMinuit;
 }

Member Function Documentation

double BCIntegrate::CalculateIntegrationVolume ( )

Calculate the integration volume

Returns: integration volume

Definition at line 248 of file BCIntegrate.cxx.

                                                {
    double integrationVolume = -1.;
 
    for(unsigned i = 0; i < fParameters.Size(); i++)
       if ( ! fParameters[i]->Fixed()) {
          if (integrationVolume<0)
             integrationVolume = 1;
          integrationVolume *= fParameters[i]->GetRangeWidth();
       }
 
    if (integrationVolume<0)
       integrationVolume = 0;
 
    return integrationVolume;
 }

bool BCIntegrate::CheckMarginalizationAvailability ( BCMarginalizationMethod type)

Check availability of integration routine for marginalization

Definition at line 558 of file BCIntegrate.cxx.

                                                                                {
   switch(type)
     {
     case BCIntegrate::kMargMonteCarlo:
       return false;
     case BCIntegrate::kMargMetropolis:
       return true;
     case BCIntegrate::kMargGrid:
       return true;
     case BCIntegrate::kMargDefault:
       return true;
     default:
       BCLog::OutError(TString::Format("BCIntegrate::CheckMarginalizationAvailability. Invalid marginalization method: %d.", type));
       break;
     }
   return false;
 }

bool BCIntegrate::CheckMarginalizationIndices	(	TH1 *	hist,
		const std::vector< unsigned > &	index
	)

Check that indices of parameters to marginalize w/r/t are correct

Definition at line 577 of file BCIntegrate.cxx.

                                                                                          {
    if (index.size()==0) {
       BCLog::OutError("BCIntegrate::Marginalize : No marginalization parameters chosen.");
       return false;
    }
 
    if (index.size() >= 4 or index.size() > fParameters.Size()) {
       BCLog::OutError("BCIntegrate::Marginalize : Too many marginalization parameters.");
       return false;
    }
 
    if ((int)index.size()<hist->GetDimension()) {
       BCLog::OutError(TString::Format("BCIntegrate::Marginalize : Too few (%d) indices supplied for histogram dimension (%d)",(int)index.size(),hist->GetDimension()));
       return false;
    }
 
    for (unsigned i=0; i<index.size(); i++) {
       // check if indices are in bounds
       if ( ! fParameters.ValidIndex(index[i])) {
          BCLog::OutError(TString::Format("BCIntegrate::Marginalize : Parameter index (%d) out of bound.",index[i]));
          return false;
       }
       // check for duplicate indices
       for (unsigned j=0; j<index.size(); j++)
          if (i!=j and index[i]==index[j]) {
             BCLog::OutError(TString::Format("BCIntegrate::Marginalize : Parameter index (%d) appears more than once",index[i]));
             return false;
          }
    }
    return true;
 }

void BCIntegrate::Copy ( const BCIntegrate & bcintegrate)

protected

Copy into object

Parameters

bcintegrate BCIntegrate object to copy values from

Definition at line 122 of file BCIntegrate.cxx.

 {
    BCEngineMCMC::Copy(other);
 
    fID                       = other.fID;
    fBestFitParameters        = other.fBestFitParameters;
    fBestFitParameterErrors   = other.fBestFitParameterErrors;
    fLogMaximum               = other.fLogMaximum;
    fMinuit                   = new TMinuit();
    fMinuitArglist[0]         = other.fMinuitArglist[0];
    fMinuitArglist[1]         = other.fMinuitArglist[1];
    fMinuitErrorFlag          = other.fMinuitErrorFlag;
    fFlagIgnorePrevOptimization = other.fFlagIgnorePrevOptimization;
    fSAT0                     = other.fSAT0;
    fSATmin                   = other.fSATmin;
    fTreeSA = 0;
    fFlagWriteSAToFile        = other.fFlagWriteSAToFile;
    fSANIterations            = other.fSANIterations;
    fSATemperature            = other.fSATemperature;
    fSALogProb                = other.fSALogProb;
    fSAx                      = other.fSAx;
    fMarginalized1D           = other.fMarginalized1D;
    fMarginalized2D           = other.fMarginalized2D;
    fOptimizationMethodCurrent= other.fOptimizationMethodCurrent;
    fOptimizationMethodUsed   = other.fOptimizationMethodUsed;
    fIntegrationMethodCurrent = other.fIntegrationMethodCurrent;
    fIntegrationMethodUsed    = other.fIntegrationMethodUsed;
    fMarginalizationMethodCurrent = other.fMarginalizationMethodCurrent;
    fMarginalizationMethodUsed= other.fMarginalizationMethodUsed;
    fSASchedule               = other.fSASchedule;
    fNIterationsMin           = other.fNIterationsMin;
    fNIterationsMax           = other.fNIterationsMax;
    fNIterationsPrecisionCheck= other.fNIterationsPrecisionCheck;
    fNIterationsOutput        = other.fNIterationsOutput;
    fNIterations              = other.fNIterations;
    fIntegral                 = other.fIntegral;
    fRelativePrecision        = other.fRelativePrecision;
    fAbsolutePrecision        = other.fAbsolutePrecision;
    fError                    = other.fError;
    fCubaIntegrationMethod    = other.fCubaIntegrationMethod;
    fCubaVegasOptions         = other.fCubaVegasOptions;
    fCubaSuaveOptions         = other.fCubaSuaveOptions;
    fCubaDivonneOptions       = other.fCubaDivonneOptions;
    fCubaCuhreOptions         = other.fCubaCuhreOptions;
    fFlagMarginalized         = other.fFlagMarginalized;
 }

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

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: An error code

Definition at line 2194 of file BCIntegrate.cxx.

 {
    BCIntegrate * local_this = static_cast<BCIntegrate *>(userdata);
 
    // scale variables
    double jacobian = 1.0;
 
    // create local parameter vector
    // important for thread safety, though not super efficient
    std::vector<double> scaled_parameters(local_this->fParameters.Size());
 
    // stay in sync with the possible lower number of parameters
    // that cuba sees due to fixing in BAT
    unsigned cubaIndex = 0;
    unsigned batIndex = 0;
    for (batIndex = 0; batIndex < local_this->fParameters.Size(); ++batIndex) {
     const BCParameter * p = local_this->fParameters[batIndex];
 
        // get the scaled parameter value
        if (p->Fixed())
           scaled_parameters[batIndex] = p->GetFixedValue();
        else {
           // convert from unit hypercube to actual parameter hyperrectangle
           scaled_parameters[batIndex] = p->GetLowerLimit() + xx[cubaIndex] * p->GetRangeWidth();
 
           // multiply range to jacobian
           jacobian *= p->GetRangeWidth();
 
           // one more parameter that cuba varies
           ++cubaIndex;
        }
    }
 
    if (cubaIndex != unsigned(*ndim))
     BCLog::OutError(Form("BCIntegrate::CubaIntegrand: mismatch between variable parameters"
                                    "in BAT (%d) and Cuba(%d)", batIndex, cubaIndex));
 
    // call function to integrate
    ff[0] = local_this->Eval(scaled_parameters);
 
    // multiply jacobian
    ff[0] *= jacobian;
 
    return 0;
 }

std::string BCIntegrate::DumpCubaIntegrationMethod ( BCIntegrate::BCCubaMethod type)

Return string with the name for a given Cuba integration type.

Parameters

type	code for the Cuba integration type

Returns: string containing the name of the Cuba integration type

Definition at line 2494 of file BCIntegrate.cxx.

 {
    switch(type) {
       case BCIntegrate::kCubaVegas:
          return "Vegas";
       case BCIntegrate::kCubaSuave:
          return "Suave";
       case BCIntegrate::kCubaDivonne:
          return "Divonne";
       case BCIntegrate::kCubaCuhre:
          return "Cuhre";
       default:
          return "Undefined";
    }
 }

std::string BCIntegrate::DumpCubaIntegrationMethod ( )

inline

Return string with the name for the currently set Cuba integration type.

Returns: string containing the name of the Cuba integration type

Definition at line 821 of file BCIntegrate.h.

822 { return DumpCubaIntegrationMethod(fCubaIntegrationMethod); }

std::string BCIntegrate::DumpCurrentIntegrationMethod ( )

inline

Return string with the name for the currently set integration type.

Returns: string containing the name of the integration type

Definition at line 767 of file BCIntegrate.h.

768 { return DumpIntegrationMethod(fIntegrationMethodCurrent); }

std::string BCIntegrate::DumpCurrentMarginalizationMethod ( )

inline

Return string with the name for the currently set marginalization type.

Returns: string containing the name of the marginalization type

Definition at line 785 of file BCIntegrate.h.

786 { return DumpMarginalizationMethod(fMarginalizationMethodCurrent); }

std::string BCIntegrate::DumpCurrentOptimizationMethod ( )

inline

Return string with the name for the currently set optimization type.

Returns: string containing the name of the optimization type

Definition at line 803 of file BCIntegrate.h.

804 { return DumpOptimizationMethod(fOptimizationMethodCurrent); }

std::string BCIntegrate::DumpIntegrationMethod ( BCIntegrate::BCIntegrationMethod type)

Return string with the name for a given integration type.

Parameters

type	code for the integration type

Returns: string containing the name of the integration type

Definition at line 2439 of file BCIntegrate.cxx.

 {
    switch(type) {
       case BCIntegrate::kIntEmpty:
          return "Empty";
       case BCIntegrate::kIntMonteCarlo:
          return "Sample Mean Monte Carlo";
       case BCIntegrate::kIntCuba:
          return "Cuba";
       case BCIntegrate::kIntGrid:
          return "Grid";
       default:
          return "Undefined";
    }
 }

std::string BCIntegrate::DumpMarginalizationMethod ( BCIntegrate::BCMarginalizationMethod type)

Return string with the name for a given marginalization type.

Parameters

type	code for the marginalization type

Returns: string containing the name of the marginalization type

Definition at line 2456 of file BCIntegrate.cxx.

 {
    switch(type) {
       case BCIntegrate::kMargEmpty:
          return "Empty";
       case BCIntegrate::kMargMonteCarlo:
          return "Sample Mean Monte Carlo";
       case BCIntegrate::kMargMetropolis:
          return "Metropolis";
       case BCIntegrate::kMargGrid:
          return "Grid";
       case BCIntegrate::kMargDefault:
          return "Default";
       default:
          return "Undefined";
    }
 }

std::string BCIntegrate::DumpOptimizationMethod ( BCIntegrate::BCOptimizationMethod type)

Return string with the name for a given optimization type.

Parameters

type	code for the optimization type

Returns: string containing the name of the optimization type

Definition at line 2475 of file BCIntegrate.cxx.

 {
    switch(type) {
       case BCIntegrate::kOptEmpty:
          return "Empty";
       case BCIntegrate::kOptSimAnn:
          return "Simulated Annealing";
       case BCIntegrate::kOptMetropolis:
          return "Metropolis MCMC";
       case BCIntegrate::kOptMinuit:
          return "Minuit";
       case BCIntegrate::kOptDefault:
          return "Default";
       default:
          return "Undefined";
    }
 }

std::string BCIntegrate::DumpUsedIntegrationMethod ( )

inline

Return string with the name for the currently set integration type.

Returns: string containing the name of the integration type

Definition at line 773 of file BCIntegrate.h.

774 { return DumpIntegrationMethod(fIntegrationMethodUsed); }

std::string BCIntegrate::DumpUsedMarginalizationMethod ( )

inline

Return string with the name for the marginalization type used.

Returns: string containing the name of the marginalization type

Definition at line 791 of file BCIntegrate.h.

792 { return DumpMarginalizationMethod(fMarginalizationMethodUsed); }

std::string BCIntegrate::DumpUsedOptimizationMethod ( )

inline

Return string with the name for the optimization type used to find the current mode.

Returns: string containing the name of the optimization type

Definition at line 809 of file BCIntegrate.h.

810 { return DumpOptimizationMethod(fOptimizationMethodUsed); }

double BCIntegrate::Eval ( const 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 211 of file BCIntegrate.cxx.

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

double BCIntegrate::EvaluatorMC	(	std::vector< double > &	sums,
		const std::vector< double > &	point,
		bool &	accepted
	)

Definition at line 535 of file BCIntegrate.cxx.

                                                                                                        {
    const double value = Eval(point);
 
    // all samples accepted in sample-mean integration
    accepted = true;
 
    // add value to sum and sum of squares
    sums[0] += value;
    sums[1] += value * value;
 
    return value;
 }

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

static

Definition at line 2115 of file BCIntegrate.cxx.

 {
    // copy parameters
    static std::vector<double> parameters;
 
    // calculate number of active + fixed parameters
    // remember: npar is just the number of _active_ parameters while
    // par is a vector of _all_ parameters
    int nparameters = ::BCIntegrateHolder::instance()->GetNParameters();
 
    // adjust size if needed
    parameters.resize(nparameters, 0.0);
 
    // copy values
    std::copy(par, par + nparameters, parameters.begin());
 
    // evaluate, for efficiency don't check if npar matches
    fval = - ::BCIntegrateHolder::instance()->LogEval(parameters);
 }

std::vector< double > BCIntegrate::FindMode ( std::vector< double > start = std::vector<double>())

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

A starting point for the mode finding can be specified for Minuit. If not specified, the previously found maximum (typically from marginalization) is used as an initial point. If that is not available, then the Minuit default will be used (center of the parameter space).

Returns: The mode found.

Note: The result may not coincide with the result of
GetBestFitParameters()

if a previous optimization found a better value.

Definition at line 1542 of file BCIntegrate.cxx.

 {
   if (fParameters.Size() < 1) {
     BCLog::OutError("FindMode : No parameters defined. Aborting.");
     return std::vector<double>();
   }
 
   if (start.empty() && fBestFitParameters.size() == fParameters.Size())
      start = fBestFitParameters;
 
   std::vector<double> mode_temp(GetNParameters());
   std::vector<double> errors_temp(GetNParameters());
   BCIntegrate::BCOptimizationMethod method_temp = fOptimizationMethodCurrent;
 
   // output
   if (!(fOptimizationMethodCurrent == BCIntegrate::kOptDefault) && !(fOptimizationMethodCurrent == BCIntegrate::kOptEmpty))
     BCLog::OutSummary(Form("Finding mode using %s", DumpCurrentOptimizationMethod().c_str()));
 
   switch (fOptimizationMethodCurrent)
     {
     case BCIntegrate::kOptEmpty:
       {
         BCLog::OutWarning("BCIntegrate::FindMode : No optimization method chosen.");
         return std::vector<double>();
       }
     case BCIntegrate::kOptSimAnn:
       {
         FindModeSA(mode_temp, errors_temp, start);
 
         break;
       }
 
     case BCIntegrate::kOptMinuit:
       {
         int printlevel = -1;
         if (BCLog::GetLogLevelScreen() <= BCLog::detail)
           printlevel = 0;
         if (BCLog::GetLogLevelScreen() <= BCLog::debug)
           printlevel = 1;
 
         BCIntegrate::FindModeMinuit(mode_temp, errors_temp, start, printlevel);
 
         break;
       }
 
     case BCIntegrate::kOptMetropolis:
       {
         FindModeMCMC(mode_temp, errors_temp);
 
         break;
       }
     case BCIntegrate::kOptDefault:
       {
         SetOptimizationMethod(BCIntegrate::kOptMinuit);
 
         return FindMode(start);
       }
 
     default:
       BCLog::OutError(Form("BCIntegrate::FindMode : Invalid mode finding method: %d", GetOptimizationMethod()));
       return std::vector<double>();
     }
 
   // calculate function at new mode
   double fcnatmode_temp = Eval(mode_temp);
 
   // check if the mode found mode is better than the previous estimate
   if ( (!fFlagIgnorePrevOptimization) && (fcnatmode_temp < fLogMaximum) ) {
     // set best fit parameters
     if (mode_temp.size() == fBestFitParameters.size())
       mode_temp      = fBestFitParameters;
     else
       BCLog::OutError("BCIntegrate::FindMode : previous mode has a different number of parameters than current.");
 
     if (errors_temp.size() == fBestFitParameterErrors.size())
       errors_temp    = fBestFitParameterErrors;
     else
       BCLog::OutError("BCIntegrate::FindMode : previous error estimate has a different number of parameters than current.");
 
     fcnatmode_temp = fLogMaximum;
     method_temp    = fOptimizationMethodUsed;
   }
 
   // set the best-fit parameters
   fBestFitParameters      = mode_temp;
   fBestFitParameterErrors = errors_temp;
   fLogMaximum             = fcnatmode_temp;
 
   // set used optimization method
   fOptimizationMethodUsed = method_temp;
 
   // return the new mode
   return mode_temp;
 
 }

std::vector< double > BCIntegrate::FindMode	(	BCIntegrate::BCOptimizationMethod	optmethod,
		std::vector< double >	start = `std::vector<double>()`
	)

Find mode using a specific method. The original method will be reset.

Parameters

optmethod	the optimization method
start	the starting point for the optimization algorithm

Returns: the mode ::vector<double> FindMode(std::vector<double> start = std::vector<double>(0));

Definition at line 1639 of file BCIntegrate.cxx.

 {
   // remember original method
   BCOptimizationMethod method_temp = fOptimizationMethodCurrent;
 
   // set method
   SetOptimizationMethod(optmethod);
 
   // run algorithm
   std::vector<double> mode = FindMode(start);
 
   // re-set original method
   SetOptimizationMethod(method_temp);
 
   // return mode
   return mode;
 }

std::vector< double > BCIntegrate::FindModeMCMC	(	std::vector< double > &	mode,
		std::vector< double > &	errors
	)

private

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

Parameters

mode	a reference to a vector holding the mode
errors	a reference to a vector holding the errors

Returns: The mode.

Note: The result may not coincide with the result of
GetBestFitParameters()

if a previous optimization found a better value.

Definition at line 2136 of file BCIntegrate.cxx.

 {
    // call PreRun
    MCMCMetropolisPreRun();
 
    // loop over all chains and find the maximum point
    double logProb = -std::numeric_limits<double>::max();
    unsigned index = 0;
    for (unsigned i = 0; i < fMCMCNChains; ++i) {
       if (MCMCGetMaximumLogProb().at(i) > logProb){
          index = i;
          logProb = MCMCGetMaximumLogProb().at(i);
       }
    }
 
    mode = MCMCGetMaximumPoint(index);
    errors.assign(fParameters.Size(),-1.);
 
    return mode;
 }

std::vector< double > BCIntegrate::FindModeMinuit	(	std::vector< double > &	mode,
		std::vector< double > &	errors,
		std::vector< double >	start = `std::vector<double>(0)`,
		int	printlevel = `1`
	)

private

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.
mode	a reference to a vector holding the mode
errors	a reference to a vector holding the errors

Returns: The mode found.

Note: The result may not coincide with the result of
GetBestFitParameters()

if a previous optimization found a better value.

Definition at line 1667 of file BCIntegrate.cxx.

 {
    if (fParameters.Size() < 1) {
       BCLog::OutError("BCIntegrate::FindModeMinuit : No parameters defined. Aborting.");
       return std::vector<double>();
    }
 
    bool have_start = true;
 
    // check start values
    if (start.size() == 0)
      have_start = false;
    else if (start.size() != fParameters.Size()) {
      have_start = false;
      BCLog::OutWarning("BCIntegrate::FindModeMinuit : Start point not valid (mismatch of dimensions), set to center.");
    }
 
    // check if point is allowed
    if (have_start) {
      for (unsigned int i = 0; i <fParameters.Size(); ++i) {
        if (!fParameters[i]->IsValid(start[i]))
          have_start = false;
        if (fParameters[i]->Fixed() && start[i] != fParameters[i]->GetFixedValue())
          have_start = false;
      }
      if (!have_start)
        BCLog::OutWarning("BCIntegrate::FindModeMinuit : Start point not valid (parameter not inside valid range), set to center.");
    }
 
    // set global this
    ::BCIntegrateHolder::instance(this);
 
    // define minuit
    if (fMinuit)
       delete fMinuit;
    fMinuit = new TMinuit(fParameters.Size());
 
    // 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 (unsigned i = 0; i < fParameters.Size(); i++) {
       double starting_point = (fParameters[i]->GetUpperLimit() + fParameters[i]->GetLowerLimit()) / 2.;
       if(have_start)
          starting_point = start[i];
       if (fParameters[i]->Fixed())
          starting_point = fParameters[i]->GetFixedValue();
       fMinuit->mnparm(i,
                       fParameters[i]->GetName().data(),
                       starting_point,
                       fParameters[i]->GetRangeWidth() / 100.,
                       fParameters[i]->GetLowerLimit(),
                       fParameters[i]->GetUpperLimit(),
                       flag);
    }
 
    for (unsigned i = 0; i < fParameters.Size(); i++)
       if (fParameters[i]->Fixed())
          fMinuit->FixParameter(i);
 
    // do mcmc minimization
    //   fMinuit->mnseek();
 
    // do minimization
    fMinuit->mnexcm("MIGRAD", fMinuitArglist, 2, flag);
 
    // improve search for local minimum
    //   fMinuit->mnimpr();
 
    // copy flag and result
    fMinuitErrorFlag = flag;
    //   std::vector<double> localMode(fParameters.Size(), 0);
    //   std::vector<double> errors(fParameters.Size(), 0);
    for (unsigned i = 0; i < fParameters.Size(); i++) {
       fMinuit->GetParameter(i, mode[i], errors[i]);
    }
 
    // delete minuit
    delete fMinuit;
    fMinuit = 0;
 
    return mode;
 }

std::vector< double > BCIntegrate::FindModeSA	(	std::vector< double > &	mode,
		std::vector< double > &	errors,
		std::vector< double >	start = `std::vector<double>(0)`
	)

private

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

Parameters

mode	a reference to a vector holding the mode
errors	a reference to a vector holding the errors
start	point in parameter space from thich the mode finding is started.

Returns: The mode.

Note: The result may not coincide with the result of
GetBestFitParameters()

if a previous optimization found a better value.

Definition at line 1777 of file BCIntegrate.cxx.

 {
    // note: if f(x) is the function to be minimized, then
    // f(x) := - LogEval(parameters)
 
    bool have_start = true;
    std::vector<double> x, y; // vectors for current state, new proposed state and best fit up to now
    double fval_x, fval_y, fval_mode; // function values at points x, y and mode (we save them rather than to re-calculate them every time)
    int t = 1; // time iterator
 
    // check start values
    if (start.size() == 0)
      have_start = false;
    else if (start.size() != fParameters.Size()) {
      have_start = false;
      BCLog::OutWarning("BCIntegrate::FindModeSA : Start point not valid (mismatch of dimensions), set to center.");
    }
 
    // check if point is allowed
    if (have_start) {
      for (unsigned int i = 0; i <fParameters.Size(); ++i) {
        if (!fParameters[i]->IsValid(start[i]))
          have_start = false;
        if (fParameters[i]->Fixed() && start[i] != fParameters[i]->GetFixedValue())
          have_start = false;
      }
      if (!have_start)
        BCLog::OutWarning("BCIntegrate::FindModeSA : Start point not valid (parameter not inside valid range), set to center.");
    }
 
    // if no starting point is given, set to center of parameter space
    if ( !have_start ) {
       start.clear();
       for (unsigned i=0; i<fParameters.Size(); i++)
          if (fParameters[i]->Fixed())
             start.push_back(fParameters[i]->GetFixedValue());
          else
             start.push_back((fParameters[i]->GetLowerLimit() + fParameters[i]->GetUpperLimit()) / 2.);
    }
 
    // set current state and best fit to starting point
    x = start;
    mode = start;
 
    // calculate function value at starting point
    fval_x = fval_mode = LogEval(x);
 
    // run while still "hot enough"
    while ( SATemperature(t) > fSATmin ) {
       // generate new state
       y = GetProposalPointSA(x, t);
 
       // check if the proposed point is inside the phase space
       // if not, reject it
       bool in_range = true;
 
       for (unsigned i = 0; i < fParameters.Size(); i++)
          if ( !fParameters[i]->IsValid(y[i]))
             in_range = false;
 
       if ( in_range ){
          // calculate function value at new point
          fval_y = LogEval(y);
 
          // is it better than the last one?
          // if so, update state and chef if it is the new best fit...
          if (fval_y >= fval_x) {
             x = y;
 
             fval_x = fval_y;
 
             if (fval_y > fval_mode) {
                mode = y;
                fval_mode = fval_y;
             }
          }
          // ...else, only accept new state w/ certain probability
          else {
             if (fRandom->Rndm() <= exp( (fval_y - fval_x) / SATemperature(t) )) {
                x = y;
                fval_x = fval_y;
             }
          }
       }
 
       // update tree variables
       fSANIterations = t;
       fSATemperature = SATemperature(t);
       fSALogProb = fval_x;
       fSAx = x;
 
       // fill tree
       if (fFlagWriteSAToFile && fTreeSA)
          fTreeSA->Fill();
 
       // increate t
       t++;
    }
 
    // calculate uncertainty
    errors.assign(fParameters.Size(),-1);
 
    return mode;
 }

double BCIntegrate::GetAbsolutePrecision ( ) const

inline

Returns: Requested absolute precision of the numerical integration

Definition at line 322 of file BCIntegrate.h.

323 { return fAbsolutePrecision; }

double BCIntegrate::GetBestFitParameter ( unsigned index) const

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

Parameters

index index of the parameter.

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

Definition at line 176 of file BCIntegrate.cxx.

 {
    if( ! fParameters.ValidIndex(index))
       return -1e+111;
 
    if(fBestFitParameters.empty()) {
       BCLog::OutError("BCIntegrate::GetBestFitParameter : Mode finding not yet run, returning center of the range.");
       return (fParameters[index]->GetUpperLimit() + fParameters[index]->GetLowerLimit() ) / 2.;
    }
 
    return fBestFitParameters[index];
 }

double BCIntegrate::GetBestFitParameterError ( unsigned index) const

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

Parameters

index index of the parameter.

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

Definition at line 190 of file BCIntegrate.cxx.

 {
    if( ! fParameters.ValidIndex(index)) {
       BCLog::OutError("BCIntegrate::GetBestFitParameterError : Parameter index out of range, returning -1.");
       return -1;
    }
 
    /*
    if(fBestFitParameterErrors.size()==0) {
       BCLog::OutError("BCIntegrate::GetBestFitParameterError : Mode finding not yet run, returning -1.");
       return -1.;
    }
 
    if(fBestFitParameterErrors[index]<0.)
       BCLog::OutWarning("BCIntegrate::GetBestFitParameterError : Parameter error not available, returning -1.");
    */
 
    return fBestFitParameterErrors[index];
 }

const std::vector<double>& BCIntegrate::GetBestFitParameterErrors ( ) const

inline

Returns the set of errors on the values of the parameters at the global mode

Definition at line 447 of file BCIntegrate.h.

448 { return fBestFitParameterErrors; }

const std::vector<double>& BCIntegrate::GetBestFitParameters ( ) const

inline

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

Returns: The best fit parameters

Definition at line 442 of file BCIntegrate.h.

443 { return fBestFitParameters; }

const BCCubaOptions::Cuhre& BCIntegrate::GetCubaCuhreOptions ( ) const

inline

Definition at line 342 of file BCIntegrate.h.

343 { return fCubaCuhreOptions; }

const BCCubaOptions::Divonne& BCIntegrate::GetCubaDivonneOptions ( ) const

inline

Definition at line 339 of file BCIntegrate.h.

340 { return fCubaDivonneOptions; }

BCCubaMethod BCIntegrate::GetCubaIntegrationMethod ( ) const

inline

Returns: Cuba Integration method

Definition at line 327 of file BCIntegrate.h.

328 { return fCubaIntegrationMethod; }

const BCCubaOptions::Suave& BCIntegrate::GetCubaSuaveOptions ( ) const

inline

Definition at line 336 of file BCIntegrate.h.

337 { return fCubaSuaveOptions; }

const BCCubaOptions::Vegas& BCIntegrate::GetCubaVegasOptions ( ) const

inline

Returns: Options used for integration with CUBA

Definition at line 333 of file BCIntegrate.h.

334 { return fCubaVegasOptions; }

double BCIntegrate::GetError ( ) const

inline

Returns: The uncertainty in the most recent Monte Carlo integration

Definition at line 396 of file BCIntegrate.h.

397 { return fError; }

double BCIntegrate::GetIntegral ( ) const

inline

Returns: The integral.

Definition at line 250 of file BCIntegrate.h.

251 { return fIntegral; }

BCIntegrate::BCIntegrationMethod BCIntegrate::GetIntegrationMethod ( ) const

inline

Returns: The current integration method

Definition at line 260 of file BCIntegrate.h.

261 { return fIntegrationMethodCurrent; }

double BCIntegrate::GetLogMaximum ( )

inline

Returns the posterior at the mode.

Returns: the posterior.

Definition at line 435 of file BCIntegrate.h.

436 { return fLogMaximum; };

BCIntegrate::BCMarginalizationMethod BCIntegrate::GetMarginalizationMethod ( ) const

inline

Returns: The current marginalization method

Definition at line 265 of file BCIntegrate.h.

266 { return fMarginalizationMethodCurrent; }

BCH1D * BCIntegrate::GetMarginalized ( const BCParameter * parameter)

Obtain the individual marginalized distributions with respect to one parameter.

Note: The most efficient method is to access by index.; Ownership of the returned heap object is conferred to the caller.

Parameters

parameter Model parameter

Returns: 1D marginalized probability

Definition at line 1200 of file BCIntegrate.cxx.

 {
    if ( !parameter)
       return 0;
    return GetMarginalized(fParameters.Index(parameter->GetName()));
 }

BCH1D* BCIntegrate::GetMarginalized ( const char * name)

inline

Obtain the individual marginalized distributions with respect to one parameter.

Note: The most efficient method is to access by index.; Ownership of the returned heap object is conferred to the caller.

Parameters

name	The parameter's name

Returns: 1D marginalized probability

Definition at line 197 of file BCIntegrate.h.

198 { return GetMarginalized(fParameters.Index(name)); }

BCH1D * BCIntegrate::GetMarginalized ( unsigned index)

Obtain the individual marginalized distributions with respect to one parameter.

Note: The most efficient method is to access by index.; Ownership of the returned heap object is conferred to the caller.

Parameters

index The parameter index

Returns: 1D marginalized probability

Definition at line 1208 of file BCIntegrate.cxx.

 {
   // check if marginalized histogram exists
   if (fMarginalized1D.size() <= index)
     return 0;
 
   // get histogram
   BCH1D * htemp = fMarginalized1D.at(index);
 
   // check if histogram exists
   if (!htemp)
     return 0;
 
   // set global mode
   if (fBestFitParameters.size() == GetNParameters())
     htemp->SetGlobalMode(fBestFitParameters.at(index));
 
   // set axis labels
    htemp->GetHistogram()->SetName(Form("hist_%i_%s", fID, fParameters[index]->GetName().data()));
    htemp->GetHistogram()->SetYTitle(Form("p(%s|data)", fParameters[index]->GetLatexName().data()));
 
    // return histogram
    return htemp;
 }

BCH2D * BCIntegrate::GetMarginalized	(	const BCParameter *	parameter1,
		const BCParameter *	parameter2
	)

Obtain the individual marginalized distributions with respect to two parameters.

Note: The most efficient method is to access by indices.; Ownership of the returned heap object is conferred to the caller.

Parameters

parameter1	First parameter
parameter2	Second parameter

Returns: 2D marginalized probability

Definition at line 1497 of file BCIntegrate.cxx.

 {
    if ( !par1 or !par2 or (par1 == par2) )
       return 0;
 
   return GetMarginalized(fParameters.Index(par1->GetName()), fParameters.Index(par2->GetName()));
 }

BCH2D* BCIntegrate::GetMarginalized	(	const char *	name1,
		const char *	name2
	)

inline

Obtain the individual marginalized distributions with respect to two parameters.

Note: The most efficient method is to access by indices.; Ownership of the returned heap object is conferred to the caller.

Parameters

name1	Name of first parameter
name2	Name of second parameter

Returns: 2D marginalized probability

Definition at line 228 of file BCIntegrate.h.

229 { return GetMarginalized(fParameters.Index(name1), fParameters.Index(name2)); }

BCH2D * BCIntegrate::GetMarginalized	(	unsigned	index1,
		unsigned	index2
	)

Obtain the individual marginalized distributions with respect to two parameters.

Note: The most efficient method is to access by indices.; Ownership of the returned heap object is conferred to the caller.

Parameters

index1	Index of first parameter
index2	Index of second parameter

Returns: 2D marginalized probability

Definition at line 1506 of file BCIntegrate.cxx.

 {
   // swap indices if necessary
   if (index1 > index2) {
     unsigned indexTemp = index1;
     index1 = index2;
     index2 = indexTemp;
   }
 
   // check if marginalized histogram exists
   if (fMarginalized2D.size() <= GetNParameters() * index1 - (index1 * index1 + 3 * index1) / 2 + index2 - 1)
     return 0;
 
   // get histogram
   BCH2D * htemp =  fMarginalized2D.at(GetNParameters() * index1 - (index1 * index1 + 3 * index1) / 2 + index2 - 1);
 
   // check if histogram exists
   if ( !htemp)
     return 0;
 
   // set global mode
   if (fBestFitParameters.size() == GetNParameters()) {
     double gmode[] = { fBestFitParameters.at(index1), fBestFitParameters.at(index2) };
     htemp->SetGlobalMode(gmode);
   }
 
   // set name
   htemp->GetHistogram()->SetName(Form("hist_%i_%s_%s", fID,
                                       fParameters[index1]->GetName().data(),
                                       fParameters[index2]->GetName().data()));
 
   // return histogram
   return htemp;
 }

TMinuit * BCIntegrate::GetMinuit ( )

Returns: Minuit used for mode finding

Definition at line 1658 of file BCIntegrate.cxx.

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

int BCIntegrate::GetMinuitErrorFlag ( ) const

inline

Returns: Error flag from Minuit run

Definition at line 405 of file BCIntegrate.h.

406 { return fMinuitErrorFlag; }

unsigned BCIntegrate::GetNIntegrationVariables ( )

Get number of variables that are varied in the integration

Returns: fNvar minus the number of fixed variables

Definition at line 239 of file BCIntegrate.cxx.

                                                {
     unsigned n = 0;
     for(unsigned i = 0; i < fParameters.Size(); ++i)
         if ( ! fParameters[i]->Fixed())
             ++n;
     return n;
 }

int BCIntegrate::GetNIterations ( ) const

inline

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

Definition at line 312 of file BCIntegrate.h.

313 { return fNIterations; }

int BCIntegrate::GetNIterationsMax ( ) const

inline

Returns: The number of maximum iterations for integration

Definition at line 297 of file BCIntegrate.h.

298 { return fNIterationsMax; }

int BCIntegrate::GetNIterationsMin ( ) const

inline

Returns: The number of minimum iterations for integration

Definition at line 292 of file BCIntegrate.h.

293 { return fNIterationsMin; }

int BCIntegrate::GetNIterationsOutput ( ) const

inline

Returns: The interval for outputting during integration

Definition at line 307 of file BCIntegrate.h.

308 { return fNIterationsOutput; }

int BCIntegrate::GetNIterationsPrecisionCheck ( ) const

inline

Returns: The interval for checking precision in integration

Definition at line 302 of file BCIntegrate.h.

303 { return fNIterationsPrecisionCheck; }

BCIntegrate::BCOptimizationMethod BCIntegrate::GetOptimizationMethod ( ) const

inline

Returns: The current optimization method

Definition at line 255 of file BCIntegrate.h.

256 { return fOptimizationMethodCurrent; }

std::vector< double > BCIntegrate::GetProposalPointSA	(	const 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 1914 of file BCIntegrate.cxx.

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

std::vector< double > BCIntegrate::GetProposalPointSABoltzmann	(	const 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 1926 of file BCIntegrate.cxx.

 {
    std::vector<double> y;
    y.clear();
    double new_val, norm;
 
    for (unsigned i = 0; i < fParameters.Size(); i++) {
       if (fParameters[i]->Fixed()) {
          y.push_back(fParameters[i]->GetFixedValue());
       }
       else {
          norm = fParameters[i]->GetRangeWidth() * SATemperature(t) / 2.;
          new_val = x[i] + norm * fRandom->Gaus();
          y.push_back(new_val);
       }
    }
 
    return y;
 }

std::vector< double > BCIntegrate::GetProposalPointSACauchy	(	const 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 1947 of file BCIntegrate.cxx.

 {
    std::vector<double> y;
    y.clear();
 
    if (fParameters.Size() == 1) {
       double cauchy, new_val, norm;
 
       if (fParameters[0]->Fixed()) {
          y.push_back(fParameters[0]->GetFixedValue());
       }
       else {
          norm = fParameters[0]->GetRangeWidth() * SATemperature(t) / 2.;
          cauchy = tan(3.14159 * (fRandom->Rndm() - 0.5));
          new_val = x[0] + norm * cauchy;
          y.push_back(new_val);
       }
    }
    else {
       // use sampling to get radial n-dim Cauchy distribution
 
       // first generate a random point uniformly distributed on a
       // fNvar-dimensional hypersphere
       y = SAHelperGetRandomPointOnHypersphere();
 
       // scale the vector by a random factor determined by the radial
       // part of the fNvar-dimensional Cauchy distribution
       double radial = SATemperature(t) * SAHelperGetRadialCauchy();
 
       // scale y by radial part and the size of dimension i in phase space
       // afterwards, move by x
       for (unsigned i = 0; i < fParameters.Size(); i++) {
          if (fParameters[i]->Fixed()) {
             y[i] = fParameters[i]->GetFixedValue(); }
          else {
             y[i] = fParameters[i]->GetRangeWidth() * y[i] * radial / 2. + x[i];
          }
       }
    }
 
    return y;
 }

std::vector< double > BCIntegrate::GetProposalPointSACustom	(	const 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 1991 of file BCIntegrate.cxx.

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

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

 {
    GetRandomVectorInParameterSpace(x);
    return Eval(x);
 }

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

Fills a vector of random numbers x[i] between fMin[i] and fMax[i] into a vector

Parameters

A	vector of doubles

Definition at line 1187 of file BCIntegrate.cxx.

 {
     GetRandomVectorUnitHypercube(x);
 
     for (unsigned i = 0; i < fParameters.Size(); i++){
        if (fParameters[i]->Fixed())
           x[i] = fParameters[i]->GetFixedValue();
        else
           x[i] = fParameters[i]->GetLowerLimit() + x[i] * fParameters[i]->GetRangeWidth();
     }
 }

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

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

Parameters

A	vector of doubles

Definition at line 1181 of file BCIntegrate.cxx.

 {
    fRandom->RndmArray(x.size(), &x.front());
 }

double BCIntegrate::GetRelativePrecision ( ) const

inline

Returns: Requested relative precision of the numerical integation

Definition at line 317 of file BCIntegrate.h.

318 { return fRelativePrecision; }

BCIntegrate::BCSASchedule BCIntegrate::GetSASchedule ( ) const

inline

Returns: The Simulated Annealing schedule

Definition at line 270 of file BCIntegrate.h.

271 { return fSASchedule; }

double BCIntegrate::GetSAT0 ( ) const

inline

Returns the Simulated Annealing starting temperature.

Definition at line 410 of file BCIntegrate.h.

411 { return fSAT0; }

double BCIntegrate::GetSATmin ( ) const

inline

Returns the Simulated Annealing threshhold temperature.

Definition at line 415 of file BCIntegrate.h.

416 { return fSATmin; }

TTree* BCIntegrate::GetSATree ( )

inline

Getter for the tree containing the Simulated Annealing chain.

Definition at line 550 of file BCIntegrate.h.

551 { return fTreeSA; }

BCH1D * BCIntegrate::GetSlice	(	const BCParameter *	parameter,
		const std::vector< double >	parameters = `std::vector<double>(0)`,
		int	bins = `0`,
		bool	flag_norm = `true`
	)

Returns a one-dimensional slice of the pdf at the point and along a specific direction.

Parameters

parameter	The model parameter along which the slice is calculated.
parameters	The point at which the other parameters are fixed.
nbins	The number of bins of the 1D-histogram.
flag_norm,:	normalize histogram to unity or not

Returns: The 1D slice.

Definition at line 1018 of file BCIntegrate.cxx.

 {
    // check if parameter exists
    if (!parameter) {
       BCLog::OutError("BCIntegrate::GetSlice : Parameter does not exist.");
       return 0;
    }
 
   // check if parameter is fixed
   if (parameter->Fixed())
     return 0;
 
    // create local copy of parameter set
    std::vector<double> parameters_temp;
    parameters_temp = parameters;
 
    // check if parameter set if defined
    if (parameters_temp.empty() && GetNParameters() == 1) {
       parameters_temp.push_back(0);
    }
    else if (parameters_temp.empty() && GetNParameters() != 1) {
       BCLog::OutError("BCIntegrate::GetSlice : No parameters defined.");
       return 0;
    }
 
    // calculate number of bins
    if (nbins <= 0)
       nbins = parameter->GetNbins();
 
    // create histogram
    TH1D * hist = new TH1D("", "", nbins, parameter->GetLowerLimit(), parameter->GetUpperLimit());
   hist->UseCurrentStyle();
 
    // set axis labels
    hist->SetXTitle(parameter->GetLatexName().data());
    if (GetNParameters() == 1)
       hist->SetYTitle(Form("p(%s|data)", parameter->GetLatexName().data()));
    else
       hist->SetYTitle(Form("p(%s|data, all other parameters fixed)", parameter->GetLatexName().data()));
    hist->SetStats(kFALSE);
 
    // fill histogram
    for (int i = 1; i <= nbins; ++i) {
       double par_temp = hist->GetBinCenter(i);
       parameters_temp[fParameters.Index(parameter->GetName())] = par_temp;
       double prob = Eval(parameters_temp);
       hist->SetBinContent(i, prob);
    }
 
    // normalize
    if (flag_norm)
       hist->Scale(1.0/hist->Integral());
 
    // set histogram
    BCH1D * hprob = new BCH1D();
    hprob->SetHistogram(hist);
 
    return hprob;
 }

BCH1D* BCIntegrate::GetSlice	(	const char *	name,
		const std::vector< double >	parameters = `std::vector<double>(0)`,
		int	nbins = `0`,
		bool	flag_norm = `true`
	)

inline

Returns a one-dimensional slice of the pdf at the point and along a specific direction.

Parameters

name	The name of the model parameter along which the slice is calculated.
parameters	The point at which the other parameters are fixed.
nbins	The number of bins of the 1D-histogram.
flag_norm,:	normalize histogram to unity or not

Returns: The 1D slice.

Definition at line 361 of file BCIntegrate.h.

362 { return GetSlice(GetParameter(name), parameters, nbins, flag_norm); }

BCH2D * BCIntegrate::GetSlice	(	const BCParameter *	parameter1,
		const BCParameter *	parameter2,
		const std::vector< double >	parameters = `std::vector<double>(0)`,
		int	bins = `0`,
		bool	flag_norm = `true`
	)

Returns a two-dimensional slice of the pdf at the point and along two specified directions.

Parameters

parameter1	The first model parameter along which the slice is calculated.
parameter2	The second model parameter along which the slice is calculated.
parameters	The point at which the other parameters are fixed.
nbins	The number of bins of the 2D-histogram.
flag_norm,:	normalize histogram to unity or not

Returns: The 2D slice.

Definition at line 1079 of file BCIntegrate.cxx.

 {
   return GetSlice(parameter1->GetName().c_str(), parameter2->GetName().c_str(), parameters, nbins, flag_norm);
 }

BCH2D * BCIntegrate::GetSlice	(	const char *	name1,
		const char *	name2,
		const std::vector< double >	parameters = `std::vector<double>(0)`,
		int	nbins = `0`,
		bool	flag_norm = `true`
	)

Returns a two-dimensional slice of the pdf at the point and along two specified directions.

Parameters

parameter1	The name of the first model parameter along which the slice is calculated.
parameter2	The name of the second model parameter along which the slice is calculated.
parameters	The point at which the other parameters are fixed.
nbins	The number of bins of the 2D-histogram.
flag_norm,:	normalize histogram to unity or not

Returns: The 2D slice.

Definition at line 1085 of file BCIntegrate.cxx.

 {
   return GetSlice(fParameters.Index(name1), fParameters.Index(name2), parameters, nbins, flag_norm);
 }

BCH2D * BCIntegrate::GetSlice	(	unsigned	index1,
		unsigned	index2,
		const std::vector< double >	parameters = `std::vector<double>(0)`,
		int	nbins = `0`,
		bool	flag_norm = `true`
	)

Returns a two-dimensional slice of the pdf at the point and along two specified directions.

Parameters

parameter1	The name of the first model parameter along which the slice is calculated.
parameter2	The name of the second model parameter along which the slice is calculated.
parameters	The point at which the other parameters are fixed.
nbins	The number of bins of the 2D-histogram.
flag_norm,:	normalize histogram to unity or not

Returns: The 2D slice.

Definition at line 1091 of file BCIntegrate.cxx.

 {
    // check if parameters exists
    if (!fParameters.ValidIndex(index1) || !fParameters.ValidIndex(index2)) {
       BCLog::OutError("BCIntegrate::GetSlice : Parameter does not exist.");
       return 0;
    }
 
    // check if parameters are fixed
    if (fParameters[index1]->Fixed() || fParameters[index2]->Fixed())
      return 0;
 
    // create local copy of parameter set
    std::vector<double> parameters_temp;
    parameters_temp = parameters;
 
    // check number of dimensions
    if (GetNParameters() < 2) {
       BCLog::OutError("BCIntegrate::GetSlice : Number of parameters need to be at least 2.");
    }
 
    // check if parameter set if defined
    if (parameters_temp.size()==0 && GetNParameters()==2) {
       parameters_temp.push_back(0);
       parameters_temp.push_back(0);
    }
    else if (parameters_temp.size()==0 && GetNParameters()>2) {
       BCLog::OutError("BCIntegrate::GetSlice : No parameters defined.");
       return 0;
    }
 
    // calculate number of bins
    const BCParameter * p1 = fParameters.Get(index1);
    const BCParameter * p2 = fParameters.Get(index2);
    unsigned nbins1, nbins2;
    if (nbins <= 0) {
       nbins1 = p1->GetNbins();
       nbins2 = p2->GetNbins();
    } else {
       nbins1 = nbins2 = nbins;
    }
 
    // create histogram
    TH2D * hist = new TH2D("", "", nbins1, p1->GetLowerLimit(), p1->GetUpperLimit(),
                           nbins2, p2->GetLowerLimit(), p2->GetUpperLimit());
    hist->UseCurrentStyle();
 
    // set axis labels
    hist->SetXTitle(Form("%s", p1->GetLatexName().data()));
    hist->SetYTitle(Form("%s", p2->GetLatexName().data()));
    hist->SetStats(kFALSE);
 
    // fill histogram
    for (unsigned int ix = 1; ix <= nbins1; ++ix) {
       for (unsigned int iy = 1; iy <= nbins2; ++iy) {
          double par_temp1 = hist->GetXaxis()->GetBinCenter(ix);
          double par_temp2 = hist->GetYaxis()->GetBinCenter(iy);
 
          parameters_temp[index1] = par_temp1;
          parameters_temp[index2] = par_temp2;
 
          double prob = Eval(parameters_temp);
          hist->SetBinContent(ix, iy, prob);
       }
    }
 
    // calculate normalization
    double norm = hist->Integral("width");
 
    // normalize
    if (flag_norm)
       hist->Scale(1.0/norm);
 
    // set histogram
    BCH2D * hprob = new BCH2D();
    hprob->SetHistogram(hist);
 
    // renormalize
    hist->Scale(norm / hist->Integral("width"));
 
    return hprob;
 }

void BCIntegrate::InitializeSATree ( )

Initialization of the tree for the Simulated Annealing

Definition at line 1759 of file BCIntegrate.cxx.

 {
   if (fTreeSA)
     delete fTreeSA;
   fTreeSA = new TTree("SATree", "SATree");
 
    fTreeSA->Branch("Iteration",      &fSANIterations,   "iteration/I");
    // todo make consistent with examples and test
    // remove because would require fixed number of parameters, and just superfluous info
 //   fTreeSA->Branch("NParameters",    &fNvar,            "parameters/I");
    fTreeSA->Branch("Temperature",    &fSATemperature,   "temperature/D");
    fTreeSA->Branch("LogProbability", &fSALogProb,       "log(probability)/D");
 
    for (unsigned i = 0; i < fParameters.Size(); ++i)
       fTreeSA->Branch(TString::Format("Parameter%i", i), &fSAx[i], TString::Format("parameter %i/D", i));
 }

void BCIntegrate::IntegralUpdaterMC	(	const std::vector< double > &	sums,
		const int &	nIterations,
		double &	integral,
		double &	absprecision
	)

static

Definition at line 549 of file BCIntegrate.cxx.

                                                                                                                                  {
    // sample mean including the volume of the parameter space
    integral = sums[2] * sums[0] / nIterations;
 
    // unbiased estimate
    absprecision = sqrt((1.0 / (nIterations - 1)) * (sums[2] * sums[2] * sums[1] / double(nIterations) - integral * integral));
 }

double BCIntegrate::Integrate ( BCIntegrationMethod intmethod)

Does the integration over the un-normalized probability.

Parameters

intmethod The integration method to used

Returns: The normalization value

Definition at line 353 of file BCIntegrate.cxx.

 {
     // remember original method
   BCIntegrationMethod method_temp = fIntegrationMethodCurrent;
 
   // set method
   SetIntegrationMethod(intmethod);
 
   // run algorithm
   double integral = Integrate();
 
   // re-set original method
   SetIntegrationMethod(method_temp);
 
   // return integral
   return integral;
 
 }

double BCIntegrate::Integrate ( )

Perform the integration

Returns: the integral

Definition at line 265 of file BCIntegrate.cxx.

 {
   // check if parameters are defined
   if (fParameters.Size() < 1) {
     BCLog::OutError("BCIntegrate::Integrate : No parameters defined. Aborting.");
     return -1.;
   }
 
   // output
   if (!(fIntegrationMethodCurrent == BCIntegrate::kIntDefault) && !(fIntegrationMethodCurrent == BCIntegrate::kIntEmpty))
     BCLog::OutSummary(Form("Integrate using %s", DumpCurrentIntegrationMethod().c_str()));
 
   switch(fIntegrationMethodCurrent)
     {
 
       // Empty
     case BCIntegrate::kIntEmpty:
       {
         BCLog::OutWarning("BCIntegrate::Integrate : No integration method chosen.");
         return 0;
       }
 
 
       // Monte Carlo Integration
     case BCIntegrate::kIntMonteCarlo:
       {
         std::vector<double> sums (2,0.0);
         sums.push_back(CalculateIntegrationVolume());
         fIntegral = Integrate(kIntMonteCarlo,
                               &BCIntegrate::GetRandomVectorInParameterSpace,
                               &BCIntegrate::EvaluatorMC,
                               &IntegralUpdaterMC,
                               sums);
 
         // set used integration method
         fIntegrationMethodUsed = BCIntegrate::kIntMonteCarlo;
 
         // return integral
         return fIntegral;
       }
 
       // CUBA library
     case BCIntegrate::kIntCuba:
       {
         fIntegral = IntegrateCuba();
 
         // set used integration method
         fIntegrationMethodUsed = BCIntegrate::kIntCuba;
 
         // return integral
         return fIntegral;
       }
 
       // CUBA library
     case BCIntegrate::kIntGrid:
       {
         fIntegral = IntegrateSlice();
 
         // set used integration method
         fIntegrationMethodUsed = BCIntegrate::kIntGrid;
 
         // return integral
         return fIntegral;
       }
 
       // default
     case BCIntegrate::kIntDefault:
       {
         if (GetNFreeParameters() <= 2)
           SetIntegrationMethod(BCIntegrate::kIntGrid);
         else
 #ifdef HAVE_CUBA_H
           SetIntegrationMethod(BCIntegrate::kIntCuba);
 #else
         SetIntegrationMethod(BCIntegrate::kIntMonteCarlo);
 #endif
         return Integrate();
       }
 
     default:
       BCLog::OutError(TString::Format("BCIntegrate::Integrate : Invalid integration method: %d", fIntegrationMethodCurrent));
       break;
     }
 
   return 0;
 }

double BCIntegrate::Integrate	(	BCIntegrationMethod	type,
		tRandomizer	randomizer,
		tEvaluator	evaluator,
		tIntegralUpdater	updater,
		std::vector< double > &	sums
	)

Does the integration over the un-normalized probability.

Parameters

type	The integration method to used (for printing out status updates by name)
randomizer	Pointer to function to choose next random point
evaluator	Pointer to function to evaluate point
updater	Pointer to function to update integral and precision
sums	Vector of doubles holding values used in integral calculation
x	An initial point to start integration routine at

Returns: The integral value

Definition at line 467 of file BCIntegrate.cxx.

 {
    LogOutputAtStartOfIntegration(type, NCubaMethods);
 
    // reset variables
    double pmax = 0.;
    double integral = 0.;
    double absprecision = 2.*fAbsolutePrecision;
    double relprecision = 2.*fRelativePrecision;
 
    std::vector<double> randx (fParameters.Size(), 0.);
 
    // how often to print out the info line to screen
    int nwrite = IntegrationOutputFrequency();
 
    // reset number of iterations
    fNIterations = 0;
    bool accepted;
 
    // iterate while number of iterations is lower than minimum number of iterations
    // or precision is not reached and the number of iterations is lower than maximum number of iterations
    while ((GetRelativePrecision() < relprecision and GetAbsolutePrecision() < absprecision and GetNIterations() < GetNIterationsMax())
          or GetNIterations() < GetNIterationsMin())
       {
 
          // get random numbers
          (this->*randomizer)(randx);
 
          // evaluate function at sampled point
          // updating sums & checking for maximum probability
 
          SetBestFitParameters(randx, (this->*evaluator)(sums,randx,accepted), pmax);
 
          // increase number of iterations
          if (accepted)
             fNIterations++;
 
          // update precisions
          if (fNIterations % fNIterationsPrecisionCheck == 0) {
             (*updater)(sums, fNIterations, integral, absprecision);
             relprecision = absprecision / integral;
          }
 
          // write status
          if (fNIterations % nwrite == 0) {
             double temp_integral;
             double temp_absprecision;
             (*updater)(sums, fNIterations, temp_integral, temp_absprecision);
             LogOutputAtIntegrationStatusUpdate(type, temp_integral, temp_absprecision, fNIterations);
          }
       }
 
    // calculate integral
    (*updater)(sums, fNIterations, integral, absprecision);
    relprecision = absprecision / integral;
 
    if (unsigned(fNIterations) >= fMCMCNIterationsMax)
       BCLog::OutWarning("BCIntegrate::Integrate: Did not converge within maximum number of iterations");
 
    // print to log
    LogOutputAtEndOfIntegration(integral, absprecision, relprecision, fNIterations);
 
    fError = absprecision;
    return integral;
 }

double BCIntegrate::IntegrateCuba ( )

inlineprivate

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

Returns: The integral

Definition at line 962 of file BCIntegrate.h.

963 { return IntegrateCuba(fCubaIntegrationMethod); }

double BCIntegrate::IntegrateCuba ( BCCubaMethod cubatype)

private

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

Parameters

Cuba	integration method to use

Returns: The integral

Definition at line 2242 of file BCIntegrate.cxx.

                                                        {
 #if HAVE_CUBA_H
    LogOutputAtStartOfIntegration(kIntCuba, cubatype);
 
    // integrand has only one component
    static const int ncomp     = 1;
 
    // don't store integration in a slot for reuse
    static const int gridno    = -1;
 
    // we don't want dumps of internal state
    static const char * statefile = "";
 
    // cuba returns info in these variables
    int fail = 0;
    int nregions = 0;
    std::vector<double> integral(ncomp, -1);
    std::vector<double> error(ncomp, -1);
    std::vector<double> prob(ncomp, -1);
 
    // reset number of iterations
    fNIterations = 0;
 
    // Cuba needs int variable
    int nIntegrationVariables = GetNIntegrationVariables();
 
    switch (cubatype) {
 
    case BCIntegrate::kCubaVegas:
       Vegas(nIntegrationVariables, ncomp,
             &BCIntegrate::CubaIntegrand, static_cast<void *>(this),
             fRelativePrecision, fAbsolutePrecision,
             fCubaVegasOptions.flags, fRandom->GetSeed(),
             fNIterationsMin, fNIterationsMax,
             fCubaVegasOptions.nstart, fCubaVegasOptions.nincrease, fCubaVegasOptions.nbatch,
             gridno, statefile,
             &fNIterations, &fail,
             &integral[0], &error[0], &prob[0]);
       break;
 
    case BCIntegrate::kCubaSuave:
       Suave(nIntegrationVariables, ncomp,
             &BCIntegrate::CubaIntegrand, static_cast<void *>(this),
             fRelativePrecision, fAbsolutePrecision,
             fCubaSuaveOptions.flags, fRandom->GetSeed(),
             fNIterationsMin, fNIterationsMax,
             fCubaSuaveOptions.nnew, fCubaSuaveOptions.flatness,
             statefile,
             &nregions, &fNIterations, &fail,
             &integral[0], &error[0], &prob[0]);
       break;
 
    case BCIntegrate::kCubaDivonne:
       if (nIntegrationVariables < 2 or nIntegrationVariables > 33)
          BCLog::OutError("BCIntegrate::IntegrateCuba(Divonne): Divonne only works in 1 < d < 34");
       else {
          // no extra info supported
          static const int ngiven = 0;
          static const int nextra = ngiven;
          Divonne(nIntegrationVariables, ncomp,
                &BCIntegrate::CubaIntegrand, static_cast<void *>(this),
                fRelativePrecision, fAbsolutePrecision,
                fCubaDivonneOptions.flags, fRandom->GetSeed(),
                fNIterationsMin, fNIterationsMax,
                fCubaDivonneOptions.key1, fCubaDivonneOptions.key2, fCubaDivonneOptions.key3,
                fCubaDivonneOptions.maxpass, fCubaDivonneOptions.border,
                fCubaDivonneOptions.maxchisq, fCubaDivonneOptions.mindeviation,
                ngiven, nIntegrationVariables /*ldxgiven*/, NULL, nextra, NULL,
                statefile,
                &nregions, &fNIterations, &fail,
                &integral[0], &error[0], &prob[0]);
       }
       break;
 
    case BCIntegrate::kCubaCuhre:
       if (nIntegrationVariables < 2)
          BCLog::OutError("BCIntegrate::IntegrateCuba(Cuhre): Cuhre(cubature) only works in d > 1");
 
       Cuhre(nIntegrationVariables, ncomp,
             &BCIntegrate::CubaIntegrand, static_cast<void *>(this),
             fRelativePrecision, fAbsolutePrecision,
             fCubaCuhreOptions.flags, fNIterationsMin, fNIterationsMax,
             fCubaCuhreOptions.key,
             statefile,
             &nregions, &fNIterations, &fail,
             &integral[0], &error[0], &prob[0]);
       break;
 
    case BCIntegrate::NCubaMethods:
    default:
       BCLog::OutError("Cuba integration method not set.");
       error[0] = -1;
       integral[0] = -1;
       break;
    }
 
    fError = error[0];
    double result = integral[0];
 
    if (fail != 0) {
       BCLog::OutWarning(" Warning, integral did not converge with the given set of parameters. ");
       BCLog::OutWarning(TString::Format(" neval    = %d", fNIterations));
       BCLog::OutWarning(TString::Format(" fail     = %d", fail));
       BCLog::OutWarning(TString::Format(" integral = %e", result));
       BCLog::OutWarning(TString::Format(" error    = %e", fError));
       BCLog::OutWarning(TString::Format(" prob     = %e", prob[0]));
 
       // handle cases in which cuba does not alter integral
       if (fNIterations == 0)
       {
          fError = -1;
          result = -1;
       }
 
    } else
       LogOutputAtEndOfIntegration(result,fError,fError/result,fNIterations);
 
    return result;
 #else
    (void) cubatype; // suppress compiler warning about unused parameters
    BCLog::OutError("IntegrateCuba: Cuba not enabled during configure");
    return -1;
 #endif
 }

double BCIntegrate::IntegrateSlice ( )

private

Integrate using the slice method

Returns: the integral;

Definition at line 2368 of file BCIntegrate.cxx.

 {
   // print to log
   LogOutputAtStartOfIntegration(fIntegrationMethodCurrent, NCubaMethods);
 
   double integral = -1;
   double absprecision  = -1;
   double relprecision  = -1;
   fError = -1;
 
   // calculate values are which the function should be evaluated
   std::vector<double> fixpoint(GetNParameters());
   for (unsigned int i = 0; i < GetNParameters(); ++i) {
     if (fParameters[i]->Fixed())
       fixpoint[i] = fParameters[i]->GetFixedValue();
     else
       fixpoint[i] = 0;
   }
 
   if (GetNFreeParameters() == 1) {
     for (unsigned int i = 0; i < GetNParameters(); ++i) {
       if (!fParameters[i]->Fixed()) {
         // calculate slice
         BCH1D* hist_temp = GetSlice(fParameters[i], fixpoint, 0, false);
 
         // calculate integral
         integral = hist_temp->GetHistogram()->Integral("width");
 
         // free memory
         delete hist_temp;
       }
     }
   }
   else if (GetNFreeParameters() == 2) {
     for (unsigned int i = 0; i < GetNParameters(); ++i) {
       for (unsigned int j = 0; j < GetNParameters(); ++j) {
         if (i >= j)
           continue;
         if (!fParameters[i]->Fixed() && !fParameters[j]->Fixed()) {
           // calculate slice
           BCH2D* hist_temp = GetSlice(GetParameter(i), GetParameter(j), fixpoint, 0, false);
 
           // calculate integral
           integral = hist_temp->GetHistogram()->Integral("width");
 
         // free memory
         delete hist_temp;
         }
       }
     }
   }
   else {
     BCLog::OutWarning("BCIntegrate::IntegrateSlice: Method only implemented for 1D and 2D problems. Return -1.");
 
     integral = -1;
   }
 
    // print to log
    LogOutputAtEndOfIntegration(integral, absprecision, relprecision, -1);
 
   return integral;
 }

unsigned BCIntegrate::IntegrationOutputFrequency ( ) const

protected

Determine frequency of output during integration

Definition at line 381 of file BCIntegrate.cxx.

 {
    if (fNIterationsOutput > 0)
       return fNIterationsOutput;
    const unsigned nwrite = fNIterationsMax / 10;
    if(nwrite < 10000)
       return 1000;
    if(nwrite < 100000)
       return 10000;
    if(nwrite < 1000000)
       return 100000;
    return 1000000;
 }

double BCIntegrate::LogEval ( const 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 218 of file BCIntegrate.cxx.

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

void BCIntegrate::LogOutputAtEndOfIntegration	(	double	integral,
		double	absprecision,
		double	relprecision,
		int	nIterations
	)

protected

Definition at line 452 of file BCIntegrate.cxx.

 {
    BCLog::OutSummary(TString::Format(" --> Result of integration:        %e +- %e", integral, absprecision));
    BCLog::OutSummary(TString::Format(" --> Obtained relative precision:  %e. ", relprecision));
    if (nIterations >= 0)
      BCLog::OutSummary(TString::Format(" --> Number of iterations:         %i", nIterations));
 }

void BCIntegrate::LogOutputAtIntegrationStatusUpdate	(	BCIntegrationMethod	type,
		double	integral,
		double	absprecision,
		int	nIterations
	)

protected

Definition at line 461 of file BCIntegrate.cxx.

 {
    BCLog::OutDetail(TString::Format("%s. Iteration %i, integral: %e +- %e.", DumpIntegrationMethod(type).c_str(), nIterations, integral, absprecision));
 }

void BCIntegrate::LogOutputAtStartOfIntegration	(	BCIntegrationMethod	type,
		BCCubaMethod	cubatype
	)

protected

Helper methods to unify output for integration methods

Parameters

type
cubatype

Definition at line 396 of file BCIntegrate.cxx.

 {
    const unsigned NVarNow = GetNIntegrationVariables();
 
    BCLog::LogLevel level = BCLog::summary;
 
    if(fParameters.Size() != NVarNow) {
 
       level = BCLog::detail;
       bool printed = false;
 
       if (type==kIntCuba)
       {
          BCLog::OutDetail(TString::Format("Running %s (%s) integration over %i dimensions out of %i.",
                DumpIntegrationMethod(type).c_str(),
                DumpCubaIntegrationMethod(cubatype).c_str(),
                NVarNow, fParameters.Size()));
          printed = true;
       }
 
       if (not printed)
          BCLog::OutDetail(TString::Format("Running %s integration over %i dimensions out of %i.",
                DumpIntegrationMethod(type).c_str(),
                NVarNow, fParameters.Size()));
 
       BCLog::OutDetail(" --> Fixed parameters:");
       for(unsigned i = 0; i < fParameters.Size(); i++)
          if(fParameters[i]->Fixed())
             BCLog::OutDetail(TString::Format("      %3i :  %g", i, fParameters[i]->GetFixedValue()));
    }
    else {
       bool printed = false;
 
       if (type==kIntCuba) {
          BCLog::OutDetail(TString::Format("Running %s (%s) integration over %i dimensions.",
                DumpIntegrationMethod(type).c_str(),
                DumpCubaIntegrationMethod(cubatype).c_str(),
                NVarNow));
          printed = true;
       }
 
       if (not printed)
          BCLog::OutDetail(TString::Format("Running %s integration over %i dimensions.",
                DumpIntegrationMethod(type).c_str(),
                NVarNow));
    }
 
    if (GetNIterationsMin() > 0 && GetNIterationsMax() > 0 ) {
      BCLog::Out(level, TString::Format(" --> Minimum number of iterations: %i", GetNIterationsMin()));
      BCLog::Out(level, TString::Format(" --> Maximum number of iterations: %i", GetNIterationsMax()));
    }
    BCLog::Out(level, TString::Format(" --> Target relative precision:    %e", GetRelativePrecision()));
    BCLog::Out(level, TString::Format(" --> Target absolute precision:    %e", GetAbsolutePrecision()));
 }

int BCIntegrate::MarginalizeAll ( )

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

Returns: Total number of marginalized distributions

Definition at line 758 of file BCIntegrate.cxx.

 {
   // check if parameters are defined
   if (fParameters.Size() < 1) {
     BCLog::OutError("BCIntegrate::MarginalizeAll : No parameters defined. Aborting.");
     return 0;
   }
 
   // check if marginalization method is defined
   if (!CheckMarginalizationAvailability(GetMarginalizationMethod())) {
     BCLog::OutError(Form("BCIntegrate::MarginalizeAll : Marginalization method not implemented \'%s\'. Aborting.", DumpCurrentMarginalizationMethod().c_str()));
     return 0;
   }
 
   // output
   if (!(fMarginalizationMethodCurrent == BCIntegrate::kMargDefault) && !(fMarginalizationMethodCurrent == BCIntegrate::kMargEmpty))
     BCLog::OutSummary(Form("Marginalize using %s", DumpCurrentMarginalizationMethod().c_str()));
 
   switch (GetMarginalizationMethod())
     {
 
       // Empty
     case BCIntegrate::kMargEmpty:
       {
         BCLog::OutWarning("BCIntegrate::MarginalizeAll : No marginalization method chosen.");
         return 0;
       }
 
       // Markov Chain Monte Carlo
     case BCIntegrate::kMargMetropolis:
       {
         // start preprocess
         MarginalizePreprocess();
 
         // run the Markov chains
         MCMCMetropolis();
 
         // start postprocess
         MarginalizePostprocess();
 
         // set pointers to marginalized distributions
         unsigned int npar = GetNParameters();
 
         fMarginalized1D.clear();
         for (unsigned int i = 0; i < npar; ++i) {
           fMarginalized1D.push_back(MCMCGetH1Marginalized(i));
         }
 
         fMarginalized2D.clear();
         if (npar > 1) {
           for (unsigned int i = 0; i < npar; ++i) {
             for (unsigned int j = 0; j < npar; ++j) {
               if (i < j)
                 fMarginalized2D.push_back(MCMCGetH2Marginalized(i, j));
             }
           }
         }
 
         // set used marginalization method
         fMarginalizationMethodUsed = BCIntegrate::kMargMetropolis;
 
         // check if mode is better than previous one
         if ( (!fFlagIgnorePrevOptimization) && (fLogMaximum < fMCMCLogMaximum) ) {
           fBestFitParameters      = fMCMCBestFitParameters;
           fBestFitParameterErrors.assign(fParameters.Size(),-1.);
           fLogMaximum             = fMCMCLogMaximum;
           fOptimizationMethodUsed = BCIntegrate::kOptMetropolis;
         }
 
       break;
       }
 
       // Sample Mean
     case BCIntegrate::kMargMonteCarlo:
       {
         return 0;
       }
 
       // Slice
     case BCIntegrate::kMargGrid:
       {
         std::vector<double> fixpoint(GetNParameters());
         for (unsigned int i = 0; i < GetNParameters(); ++i) {
           if (fParameters[i]->Fixed())
             fixpoint[i] = fParameters[i]->GetFixedValue();
           else
             fixpoint[i] = 0;
         }
 
         // start preprocess
         MarginalizePreprocess();
 
         // store highest probability
         double probmax = -std::numeric_limits<double>::infinity();
         std::vector<double> bestfit_parameters(fParameters.Size(), -std::numeric_limits<double>::infinity());
         // correct fixed parameter values
         // the other should have been determined above
         for (unsigned i = 0; i < fParameters.Size(); ++i)
            if (fParameters[i]->Fixed())
               bestfit_parameters[i] = fParameters[i]->GetFixedValue();
 
         fMarginalized1D.clear();
         if (GetNFreeParameters() == 1) {
           for (unsigned int i = 0; i < GetNParameters(); ++i) {
             if (!fParameters[i]->Fixed()) {
               // calculate slice
               BCH1D* hist_temp = GetSlice(fParameters[i], fixpoint, 0);
               fMarginalized1D.push_back(hist_temp);
 
               // remember best fit before it is normalized and the value is wrong
               const int index = hist_temp->GetHistogram()->GetMaximumBin();
               probmax = hist_temp->GetHistogram()->GetBinContent(index);
               bestfit_parameters.at(i) = hist_temp->GetHistogram()->GetBinCenter(index);
 
               // normalize to unity
               hist_temp->GetHistogram()->Scale(1./hist_temp->GetHistogram()->Integral());
             }
             else
               fMarginalized1D.push_back(0);
           }
 
           fMarginalized2D.clear();
           for (unsigned int i = 0; i < GetNParameters(); ++i) {
             for (unsigned int j = 0; j < GetNParameters(); ++j) {
               if (i >= j)
                 continue;
               else
                 fMarginalized2D.push_back(0);
             }
           }
         }
         else if (GetNFreeParameters() == 2) {
           // create a 2D histogram
 
           BCH2D* hist2d_temp = 0;
 
           fMarginalized2D.clear();
           for (unsigned int i = 0; i < GetNParameters(); ++i) {
             for (unsigned int j = 0; j < GetNParameters(); ++j) {
               if (i >= j)
                 continue;
               if (!fParameters[i]->Fixed() && !fParameters[j]->Fixed()) {
                 // calculate slice
                 hist2d_temp = GetSlice(GetParameter(i), GetParameter(j), fixpoint, 0);
                 fMarginalized2D.push_back(hist2d_temp);
 
                 // remember best fit before it is normalized and the value is wrong
                 int index1, index2, useless_index;
                 hist2d_temp->GetHistogram()->GetMaximumBin(index1, index2, useless_index);
                 probmax = hist2d_temp->GetHistogram()->GetBinContent(index1, index2, useless_index);
                 bestfit_parameters.at(i) = hist2d_temp->GetHistogram()->GetXaxis()->GetBinCenter(index1);
                 bestfit_parameters.at(j) = hist2d_temp->GetHistogram()->GetYaxis()->GetBinCenter(index2);
 
                 // normalize to unity
                 hist2d_temp->GetHistogram()->Scale(1./hist2d_temp->GetHistogram()->Integral());
               }
               else
                 fMarginalized2D.push_back(0);
             }
           }
 
           // marginalize by projecting
           fMarginalized1D.clear();
           for (unsigned int i = 0; i < GetNParameters(); ++i)
             fMarginalized1D.push_back(0);
           for (unsigned int i = 0; i < GetNParameters(); ++i) {
             for (unsigned int j = 0; j < GetNParameters(); ++j) {
               if (i >= j)
                 continue;
               if (!fParameters[i]->Fixed() && !fParameters[j]->Fixed()) {
                 BCH1D* hist_x = new BCH1D( hist2d_temp->GetHistogram()->ProjectionX() );
                 hist_x->GetHistogram()->UseCurrentStyle(); // reset style because histogram is created by a projection
                 hist_x->GetHistogram()->SetStats(kFALSE);
                 fMarginalized1D[i] = hist_x;
                 BCH1D* hist_y = new BCH1D( hist2d_temp->GetHistogram()->ProjectionY() );
                 hist_y->GetHistogram()->UseCurrentStyle(); // reset style because histogram is created by a projection
                 hist_y->GetHistogram()->SetStats(kFALSE);
                 fMarginalized1D[j] = hist_y;
               }
             }
           }
         }
         else {
           BCLog::OutWarning("BCIntegrate::MarginalizeAll : Option works for 1D and 2D problems.");
           return 0;
         }
 
         // save if improved the log posterior
 
         if (fBestFitParameters.empty() || probmax > fMCMCLogMaximum) {
            fLogMaximum = probmax;
            fBestFitParameters = bestfit_parameters;
            fBestFitParameterErrors.assign(fBestFitParameters.size(), std::numeric_limits<double>::infinity());
            // todo can't set this one yet
            //           fOptimizationMethodUsed =
         }
 
         BCLog::OutDetail(" --> Global mode from Grid:");
         BCLog::OutDebug(Form(" --> Posterior value: %g", probmax));
         int ndigits = (int) log10(fParameters.Size());
         for (unsigned i = 0; i < fParameters.Size(); ++i)
             BCLog::OutDetail(TString::Format(" -->      parameter %*i:   %.4g", ndigits+1, i, bestfit_parameters[i]));
 
         // start postprocess
         MarginalizePostprocess();
 
         // set used marginalization method
         fMarginalizationMethodUsed = BCIntegrate::kMargGrid;
 
       break;
       }
 
       // default
     case BCIntegrate::kMargDefault:
       {
         if ( GetNFreeParameters() <= 2 ) {
           SetMarginalizationMethod(BCIntegrate::kMargGrid);
         }
         else {
           SetMarginalizationMethod(BCIntegrate::kMargMetropolis);
         }
 
         // call again
         return MarginalizeAll();
 
       break;
       }
 
     default:
       BCLog::OutError(TString::Format("BCIntegrate::MarginalizeAll : Invalid marginalization method: %d", GetMarginalizationMethod()));
       return 0;
       break;
     }
 
   // set flag
   fFlagMarginalized = true;
 
   return 1;
 }

int BCIntegrate::MarginalizeAll ( BCIntegrate::BCMarginalizationMethod margmethod)

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

Parameters

margmethod the marginalization method.

Returns: Total number of marginalized distributions

Definition at line 999 of file BCIntegrate.cxx.

 {
   // remember original method
   BCMarginalizationMethod method_temp = fMarginalizationMethodCurrent;
 
   // set method
   SetMarginalizationMethod(margmethod);
 
   // run algorithm
   int result = MarginalizeAll();
 
   // re-set original method
   SetMarginalizationMethod(method_temp);
 
   // return result
   return result;
 }

virtual void BCIntegrate::MarginalizePostprocess ( )

inlinevirtual

Method executed after marginalization. User's code should be provided via overloading in the derived class

Reimplemented in BCFitter.

Definition at line 647 of file BCIntegrate.h.

648 {};

virtual void BCIntegrate::MarginalizePreprocess ( )

inlinevirtual

Method executed for before marginalization. User's code should be provided via overloading in the derived class

Reimplemented in BCFitter.

Definition at line 641 of file BCIntegrate.h.

642 {};

double BCIntegrate::Normalize ( )

inline

Performs integration.

Definition at line 578 of file BCIntegrate.h.

579 { return Integrate(); };

BCIntegrate & BCIntegrate::operator= ( const BCIntegrate & bcintegrate)

Defaut assignment operator

Definition at line 115 of file BCIntegrate.cxx.

 {
    Copy(other);
    return *this;
 }

int BCIntegrate::PrintAllMarginalized	(	const char *	file,
		std::string	options1d = `"BTsiB3CS1D0pdf0Lmeanmode"`,
		std::string	options2d = `"BTfB3CS1meangmode"`,
		unsigned int	hdiv = `1`,
		unsigned int	ndiv = `1`
	)

Definition at line 1346 of file BCIntegrate.cxx.

 {
    if (fMarginalized1D.size() == 0 and fMarginalized2D.size() == 0) {
       BCLog::OutError("BCIntegrate::PrintAllMarginalized : Marginalized distributions not available.");
       return 0;
    }
 
    // count all valid (non NULL) histograms
    int nvalid1D = 0;
    int nvalid2D = 0;
 
    for (unsigned i = 0 ; i < fMarginalized1D.size() ; ++i) {
      if (BCH1D * h = fMarginalized1D[i])
        if (h->GetHistogram())
          nvalid1D++;
    }
 
    for (unsigned i = 0 ; i < fMarginalized2D.size() ; ++i) {
      if (BCH2D * h = fMarginalized2D[i])
        if (h->GetHistogram())
          nvalid2D++;
    }
 
    std::string filename(file);
 
    // check if file extension does not exist or is not pdf or ps
    if ( (filename.find_last_of(".") == std::string::npos) or
          ((filename.substr(filename.find_last_of(".")+1) != "pdf") and
                (filename.substr(filename.find_last_of(".")+1) != "ps"))) {
       // make it a PDF file
       filename += ".pdf";
    }
 
    int c_width  = gStyle->GetCanvasDefW(); // default canvas width
    int c_height = gStyle->GetCanvasDefH(); // default canvas height
 
    if (hdiv > vdiv) {
       if (hdiv > 3) {
          c_width = 1000;
          c_height = 700;
       }
       else {
          c_width = 800;
          c_height = 600;
       }
    }
    else if (hdiv < vdiv) {
       if (hdiv > 3) {
          c_height = 1000;
          c_width = 700;
       }
       else {
          c_height = 800;
          c_width = 600;
       }
    }
 
    const unsigned nplots = nvalid1D + nvalid2D;
 
    // give out warning if too many plots
    BCLog::OutSummary(Form("Printing all marginalized distributions (%d x 1D + %d x 2D = %d) into file %s",
                           nvalid1D, nvalid2D, nplots, filename.c_str()));
    if (nplots > 100)
       BCLog::OutDetail("This can take a while...");
 
    // setup the canvas and file
    TCanvas c("c", "canvas", c_width, c_height);
    c.Divide(hdiv, vdiv);
 
    // for clean up later
    std::vector<BCH1D *> h1;
 
    // count plots
    unsigned n = 0;
    for (unsigned i = 0; i < fParameters.Size(); ++i) {
       BCH1D * h = GetMarginalized(i);
       h1.push_back(h);
 
       // check if histogram exists
       if ( !h)
          continue;
 
       // if current page is full, switch to new page
       if (n != 0 && n % (hdiv * vdiv) == 0) {
          if ( n <= (hdiv * vdiv)) {
             c.Print(std::string( filename + "(").c_str());
          }
          else {
             c.Print(filename.c_str());
          }
       }
 
       // go to next pad
       c.cd(n % (hdiv * vdiv) + 1);
 
       h->Draw(options1d);
 
       if (++n % 100 == 0)
          BCLog::OutDetail(Form(" --> %d plots done", n));
    }
 
    // for clean up later
    std::vector<BCH2D *> h2;
 
    // check how many 2D plots are actually drawn, despite no histogram filling or delta prior
    unsigned k = 0;
    for (unsigned i = 0; i < fParameters.Size() - 1; ++i) {
       for (unsigned j = i + 1; j < fParameters.Size(); ++j) {
          // check if histogram exists, or skip if one par has a delta prior
          h2.push_back(GetMarginalized(i, j));
          if ( !h2.back())
             continue;
 
          // get corresponding parameters
          BCParameter * a = GetParameter(i);
          BCParameter * b = GetParameter(j);
 
          // if current page is full, switch to new page, but only if there is data to plot
          if ((k != 0 && k % (hdiv * vdiv) == 0) || k == 0) {
             c.Print(filename.c_str());
          }
 
          // go to next pad
          c.cd(k % (hdiv * vdiv) + 1);
 
          const double meana = (a->GetLowerLimit() + a->GetUpperLimit()) / 2.;
          if (a->GetRangeWidth() <= 1e-7 * meana)
             continue;
 
          const double meanb = (b->GetLowerLimit() + b->GetUpperLimit()) / 2.;
          if (b->GetRangeWidth() <= 1e-7 * meanb)
             continue;
 
          h2.back()->Draw(options2d);
          k++;
 
          if ((n + k) % 100 == 0)
             BCLog::OutDetail(Form(" --> %d plots done", n + k));
       }
    }
 
    if ((n + k) > 100 && (n + k) % 100 != 0)
       BCLog::OutDetail(Form(" --> %d plots done", n + k));
 
    c.Print(std::string( filename + ")").c_str());
 
    // return total number of drawn histograms
    return n + k;
 }

int BCIntegrate::PrintAllMarginalized1D ( const char * filebase)

Definition at line 1293 of file BCIntegrate.cxx.

 {
    if (fMCMCH1Marginalized.size() == 0) {
       BCLog::OutError("BCIntegrate::PrintAllMarginalized : Marginalized distributions not available.");
       return 0;
    }
 
    int n = GetNParameters();
    for (int i = 0; i < n; i++) {
       const BCParameter * a = GetParameter(i);
       if (GetMarginalized(a))
          GetMarginalized(a)->Print(Form("%s_1D_%s.pdf", filebase, a->GetName().data()));
    }
 
    return n;
 }

int BCIntegrate::PrintAllMarginalized2D ( const char * filebase)

Definition at line 1311 of file BCIntegrate.cxx.

 {
    if (fMCMCH2Marginalized.size() == 0) {
       BCLog::OutError("BCIntegrate::PrintAllMarginalized : Marginalized distributions not available.");
       return 0;
    }
 
    int k = 0;
    int n = GetNParameters();
    for (int i = 0; i < n - 1; i++) {
       for (int j = i + 1; j < n; j++) {
          const BCParameter * a = GetParameter(i);
          const BCParameter * b = GetParameter(j);
 
          double meana = (a->GetLowerLimit() + a->GetUpperLimit()) / 2.;
          double deltaa = (a->GetUpperLimit() - a->GetLowerLimit());
          if (deltaa <= 1e-7 * meana)
             continue;
 
          double meanb = (b->GetLowerLimit() + b->GetUpperLimit()) / 2.;
          double deltab = (b->GetUpperLimit() - b->GetLowerLimit());
          if (deltab <= 1e-7 * meanb)
             continue;
 
          if (GetMarginalized(a, b))
             GetMarginalized(a, b)->Print(
                   Form("%s_2D_%s_%s.pdf",filebase, a->GetName().data(), b->GetName().data()));
          k++;
       }
    }
 
    return k;
 }

int BCIntegrate::ReadMarginalizedFromFile ( const char * file)

Read

Definition at line 1234 of file BCIntegrate.cxx.

 {
    TFile * froot = TFile::Open(file);
    if (!froot->IsOpen()) {
       BCLog::OutError(Form("BCIntegrate::ReadMarginalizedFromFile : Couldn't open file %s.", file));
       return 0;
    }
 
    // We reset the MCMCEngine here for the moment.
    // In the future maybe we only want to do this if the engine
    // wans't initialized at all or when there were some changes
    // in the model.
    // But maybe we want reset everything since we're overwriting
    // the marginalized distributions anyway.
    MCMCInitialize();
 
    int k = 0;
    int n = GetNParameters();
    for (int i = 0; i < n; i++) {
       const BCParameter * a = GetParameter(i);
       TKey * key = froot->GetKey(TString::Format("hist_%i_%s", fID, a->GetName().data()));
       if (key) {
          TH1D * h1 = (TH1D*) key->ReadObjectAny(TH1D::Class());
          h1->SetDirectory(0);
          if (SetMarginalized(i, h1))
             k++;
       }
       else
          BCLog::OutWarning(
                Form("BCIntegrate::ReadMarginalizedFromFile : Couldn't read histogram \"hist_%i_%s\" from file %s.",
                      fID, a->GetName().data(), file));
    }
 
    for (int i = 0; i < n - 1; i++) {
       for (int j = i + 1; j < n; j++) {
          const BCParameter * a = GetParameter(i);
          const BCParameter * b = GetParameter(j);
          TKey * key = froot->GetKey(
                TString::Format("hist_%i_%s_%s",fID, a->GetName().data(),b->GetName().data()));
          if (key) {
             TH2D * h2 = (TH2D*) key->ReadObjectAny(TH2D::Class());
             h2->SetDirectory(0);
             if (SetMarginalized(i, j, h2))
                k++;
          }
          else
             BCLog::OutWarning(
                   Form("BCIntegrate::ReadMarginalizedFromFile : Couldn't read histogram \"hist_%i_%s_%s\" from file %s.",
                         fID, a->GetName().data(), b->GetName().data(), file));
       }
    }
 
    froot->Close();
 
    return k;
 }

void BCIntegrate::ResetResults ( )

virtual

Reset all information on the best fit parameters.

Reimplemented from BCEngineMCMC.

Definition at line 225 of file BCIntegrate.cxx.

 {
    BCEngineMCMC::ResetResults();
 
    fBestFitParameterErrors.clear();
    fBestFitParameters.clear();
    fLogMaximum = -std::numeric_limits<double>::max();
 
    // remove marginalized histograms
    // set marginalization flag
    fFlagMarginalized = false;
 }

double BCIntegrate::SAHelperGetRadialCauchy ( )

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

Definition at line 2042 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 unsigned map_dimension = 0;
 
    // is the lookup-table already initialized? if not, do it!
    if (!initialized or map_dimension != fParameters.Size()) {
       double init_theta;
       double beta = SAHelperSinusToNIntegral(fParameters.Size() - 1, 1.57079632679);
 
       for (int i = 0; i <= 10000; i++) {
          init_theta = 3.14159265 * (double)i / 5000.;
          map_theta[i] = init_theta;
 
          map_u[i] = SAHelperSinusToNIntegral(fParameters.Size() - 1, init_theta) / beta;
       }
 
       map_dimension = fParameters.Size();
       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 1998 of file BCIntegrate.cxx.

 {
    std::vector<double> rand_point(fParameters.Size());
 
    // 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 (fParameters.Size() == 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 (fParameters.Size() == 3) {
       fRandom->Sphere(rand_point[0], rand_point[1], rand_point[2], 1.0);
    }
    else {
       double s = 0.,
          gauss_num;
 
       for (unsigned i = 0; i < fParameters.Size(); i++) {
          gauss_num = fRandom->Gaus();
          rand_point[i] = gauss_num;
          s += gauss_num * gauss_num;
       }
       s = sqrt(s);
 
       for (unsigned i = 0; i < fParameters.Size(); 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 2098 of file BCIntegrate.cxx.

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

void BCIntegrate::SAInitialize ( )

Initializes the Simulated Annealing algorithm (for details see manual)

Definition at line 2432 of file BCIntegrate.cxx.

 {
    fSAx.clear();
    fSAx.assign(fParameters.Size(), 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 1883 of file BCIntegrate.cxx.

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

double BCIntegrate::SATemperatureBoltzmann ( double t)

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

Parameters

t	iterator for lowering the temperature over time.

Definition at line 1895 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 1901 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 1907 of file BCIntegrate.cxx.

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

void BCIntegrate::SetAbsolutePrecision ( double absprecision)

inline

Set absolute precision of the numerical integation

Definition at line 513 of file BCIntegrate.h.

514 { fAbsolutePrecision = absprecision; }

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

inline

Set best fit parameters values

Definition at line 826 of file BCIntegrate.h.

827 { fBestFitParameters = x; }

void BCIntegrate::SetBestFitParameters	(	const std::vector< double > &	x,
		const double &	new_value,
		double &	old_value
	)

Set best fit parameters if best fit

Parameters

new_value	is the value of the function at x
old_value	is the old best fit value, updated to new_value, if it is larger

Definition at line 373 of file BCIntegrate.cxx.

                                                                                                              {
    if (new_value < old_value)
       return;
    old_value = new_value;
    SetBestFitParameters(x);
 }

void BCIntegrate::SetCubaIntegrationMethod ( BCIntegrate::BCCubaMethod type)

Set Cuba integration method

Definition at line 2173 of file BCIntegrate.cxx.

 {
 #ifdef HAVE_CUBA_H
    switch(type) {
       case BCIntegrate::kCubaVegas:
       case BCIntegrate::kCubaSuave:
       case BCIntegrate::kCubaDivonne:
       case BCIntegrate::kCubaCuhre:
          fCubaIntegrationMethod = type;
          return;
       default:
          BCLog::OutError(TString::Format("Integration method of type %d is not defined for Cuba",type));
          return;
    }
 #else
    (void) type; // suppress compiler warning about unused parameters
    BCLog::OutError("SetCubaIntegrationMethod: Cuba not enabled during configure");
 #endif
 }

void BCIntegrate::SetCubaOptions ( const BCCubaOptions::Vegas & options)

inline

Set options for individual cuba methods

Definition at line 523 of file BCIntegrate.h.

524 { fCubaVegasOptions = options; }

void BCIntegrate::SetCubaOptions ( const BCCubaOptions::Suave & options)

inline

Definition at line 526 of file BCIntegrate.h.

527 { fCubaSuaveOptions = options; }

void BCIntegrate::SetCubaOptions ( const BCCubaOptions::Divonne & options)

inline

Definition at line 529 of file BCIntegrate.h.

530 { fCubaDivonneOptions = options; }

void BCIntegrate::SetCubaOptions ( const BCCubaOptions::Cuhre & options)

inline

Definition at line 532 of file BCIntegrate.h.

533 { fCubaCuhreOptions = options; }

void BCIntegrate::SetFlagIgnorePrevOptimization ( bool flag)

inline

Flag whether or not to ignore result of previous mode finding

Definition at line 463 of file BCIntegrate.h.

464 { fFlagIgnorePrevOptimization = flag; }

void BCIntegrate::SetFlagWriteSAToFile ( bool flag)

inline

Definition at line 545 of file BCIntegrate.h.

546 { fFlagWriteSAToFile = flag; }

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

Parameters

method The current integration method

Definition at line 2158 of file BCIntegrate.cxx.

 {
    if (method >= BCIntegrate::NIntMethods) {
       BCLog::OutError(Form("BCIntegrate::SetIntegrationMethod: Invalid method '%d' ", method));
       return;
    }
    if (method == BCIntegrate::kIntCuba) {
 #ifndef HAVE_CUBA_H
       BCLog::OutError("BCIntegrate::SetIntegrationMethod: Cuba not enabled during configure");
 #endif
    }
    fIntegrationMethodCurrent = method;
 }

void BCIntegrate::SetMarginalizationMethod ( BCIntegrate::BCMarginalizationMethod method)

inline

Parameters

method The current marginalization method

Definition at line 477 of file BCIntegrate.h.

478 { fMarginalizationMethodCurrent = method; }

void BCIntegrate::SetMinuitArlist ( double * arglist)

inline

pointer to list of doubles to be passed as arguments to Minuit

Definition at line 457 of file BCIntegrate.h.

458 { fMinuitArglist[0] = arglist[0];

459 fMinuitArglist[1] = arglist[1]; }

void BCIntegrate::SetNIterationsMax ( int niterations)

inline

Parameters

niterations The maximum number of iterations for integration

Definition at line 492 of file BCIntegrate.h.

493 { fNIterationsMax = niterations; }

void BCIntegrate::SetNIterationsMin ( int niterations)

inline

Parameters

niterations The maximum number of iterations for integration

Definition at line 487 of file BCIntegrate.h.

488 { fNIterationsMin = niterations; }

void BCIntegrate::SetNIterationsOutput ( int niterations)

inline

Parameters

niterations interval for outputting during integration. If negative, frequency is autogenerated.

Definition at line 502 of file BCIntegrate.h.

503 { fNIterationsOutput = niterations; }

void BCIntegrate::SetNIterationsPrecisionCheck ( int niterations)

inline

Parameters

niterations interval for checking precision in integration routines

Definition at line 497 of file BCIntegrate.h.

498 { fNIterationsPrecisionCheck = niterations; }

void BCIntegrate::SetOptimizationMethod ( BCIntegrate::BCOptimizationMethod method)

inline

Parameters

method The current optimization method

Definition at line 468 of file BCIntegrate.h.

469 { fOptimizationMethodCurrent = method; }

void BCIntegrate::SetRelativePrecision ( double relprecision)

inline

Parameters

relprecision The relative precision envisioned for Monte Carlo integration

Definition at line 508 of file BCIntegrate.h.

509 { fRelativePrecision = relprecision; }

void BCIntegrate::SetSASchedule ( BCIntegrate::BCSASchedule schedule)

inline

Parameters

method The Simulated Annealing schedule

Definition at line 482 of file BCIntegrate.h.

483 { fSASchedule = schedule; }

void BCIntegrate::SetSAT0 ( double T0)

inline

Parameters

T0	new value for Simulated Annealing starting temperature.

Definition at line 537 of file BCIntegrate.h.

538 { fSAT0 = T0; }

void BCIntegrate::SetSATmin ( double Tmin)

inline

Parameters

Tmin	new value for Simulated Annealing threshold temperature.

Definition at line 542 of file BCIntegrate.h.

543 { fSATmin = Tmin; }

Member Data Documentation

double BCIntegrate::fAbsolutePrecision

private

Requested relative precision of the integration

Definition at line 1045 of file BCIntegrate.h.

std::vector<double> BCIntegrate::fBestFitParameterErrors

private

A vector of estimates on the uncertainties

Definition at line 1031 of file BCIntegrate.h.

std::vector<double> BCIntegrate::fBestFitParameters

private

A vector of best fit parameters found by MCMC

Definition at line 1027 of file BCIntegrate.h.

BCCubaOptions::Cuhre BCIntegrate::fCubaCuhreOptions

private

Definition at line 1056 of file BCIntegrate.h.

BCCubaOptions::Divonne BCIntegrate::fCubaDivonneOptions

private

Definition at line 1055 of file BCIntegrate.h.

BCCubaMethod BCIntegrate::fCubaIntegrationMethod

private

Cuba integration method

Definition at line 1052 of file BCIntegrate.h.

BCCubaOptions::Suave BCIntegrate::fCubaSuaveOptions

private

Definition at line 1054 of file BCIntegrate.h.

BCCubaOptions::Vegas BCIntegrate::fCubaVegasOptions

private

Definition at line 1053 of file BCIntegrate.h.

double BCIntegrate::fError

private

The uncertainty in the most recent Monte Carlo integration

Definition at line 1049 of file BCIntegrate.h.

bool BCIntegrate::fFlagIgnorePrevOptimization

protected

Flag for ignoring older results of optimization

Definition at line 870 of file BCIntegrate.h.

bool BCIntegrate::fFlagMarginalized

protected

flag indicating if the model was marginalized

Definition at line 922 of file BCIntegrate.h.

bool BCIntegrate::fFlagWriteSAToFile

protected

Flag deciding whether to write SA to file or not.

Definition at line 886 of file BCIntegrate.h.

int BCIntegrate::fID

protected

An identification number in case several models exist .

Definition at line 859 of file BCIntegrate.h.

double BCIntegrate::fIntegral

private

The integral.

Definition at line 1039 of file BCIntegrate.h.

BCIntegrate::BCIntegrationMethod BCIntegrate::fIntegrationMethodCurrent

private

Current integration method

Definition at line 987 of file BCIntegrate.h.

BCIntegrate::BCIntegrationMethod BCIntegrate::fIntegrationMethodUsed

private

Integration method used for the current results

Definition at line 991 of file BCIntegrate.h.

double BCIntegrate::fLogMaximum

private

The function value at the mode on the log scale

Definition at line 1035 of file BCIntegrate.h.

BCIntegrate::BCMarginalizationMethod BCIntegrate::fMarginalizationMethodCurrent

private

Current marginalization method

Definition at line 995 of file BCIntegrate.h.

BCIntegrate::BCMarginalizationMethod BCIntegrate::fMarginalizationMethodUsed

private

Marginalization method used for the current results

Definition at line 999 of file BCIntegrate.h.

std::vector<BCH1D*> BCIntegrate::fMarginalized1D

protected

Set of marginalized distributions.

Definition at line 895 of file BCIntegrate.h.

std::vector<BCH2D*> BCIntegrate::fMarginalized2D

protected

Set of marginalized distributions.

Definition at line 899 of file BCIntegrate.h.

TMinuit* BCIntegrate::fMinuit

protected

Minuit

Definition at line 863 of file BCIntegrate.h.

double BCIntegrate::fMinuitArglist[2]

protected

Definition at line 865 of file BCIntegrate.h.

int BCIntegrate::fMinuitErrorFlag

protected

Definition at line 866 of file BCIntegrate.h.

int BCIntegrate::fNIterations

private

Number of iterations in the most recent Monte Carlo integration

Definition at line 1023 of file BCIntegrate.h.

unsigned BCIntegrate::fNIterationsMax

private

Maximum number of iterations

Definition at line 1011 of file BCIntegrate.h.

unsigned BCIntegrate::fNIterationsMin

private

Maximum number of iterations

Definition at line 1007 of file BCIntegrate.h.

unsigned BCIntegrate::fNIterationsOutput

private

Output frequency during integration

Definition at line 1019 of file BCIntegrate.h.

unsigned BCIntegrate::fNIterationsPrecisionCheck

private

Maximum number of iterations

Definition at line 1015 of file BCIntegrate.h.

BCIntegrate::BCOptimizationMethod BCIntegrate::fOptimizationMethodCurrent

private

Current mode finding method

Definition at line 978 of file BCIntegrate.h.

BCIntegrate::BCOptimizationMethod BCIntegrate::fOptimizationMethodUsed

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

double BCIntegrate::fRelativePrecision

private

Requested relative precision of the integration

Definition at line 1042 of file BCIntegrate.h.

double BCIntegrate::fSALogProb

protected

Definition at line 890 of file BCIntegrate.h.

int BCIntegrate::fSANIterations

protected

Definition at line 888 of file BCIntegrate.h.

BCIntegrate::BCSASchedule BCIntegrate::fSASchedule

private

Current Simulated Annealing schedule

Definition at line 1003 of file BCIntegrate.h.

double BCIntegrate::fSAT0

protected

Starting temperature for Simulated Annealing

Definition at line 874 of file BCIntegrate.h.

double BCIntegrate::fSATemperature

protected

Definition at line 889 of file BCIntegrate.h.

double BCIntegrate::fSATmin

protected

Minimal/Threshold temperature for Simulated Annealing

Definition at line 878 of file BCIntegrate.h.

std::vector<double> BCIntegrate::fSAx

protected

Definition at line 891 of file BCIntegrate.h.

TTree* BCIntegrate::fTreeSA

protected

Tree for the Simulated Annealing

Definition at line 882 of file BCIntegrate.h.

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

Public Types

Public Member Functions

Protected Member Functions

Protected Attributes

Private Member Functions

Private Attributes

Member functions (miscellaneous methods)

Detailed Description

Member Typedef Documentation

Member Enumeration Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation