a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
Public Member Functions | Private Member Functions | Private Attributes | List of all members
MonteCarloEngine Class Reference
MonteCarlo

An engine class for Monte Carlo. More...

#include <MonteCarloEngine.h>

+ Inheritance diagram for MonteCarloEngine:

Detailed Description

An engine class for Monte Carlo.

Author
HEPfit Collaboration
Copyright
GNU General Public License

The MonteCarlo engine is has overloaded from BCEngineMCMC class in the BAT libraries and user defined fumctions to facilitate the Monte Carlo run.

Definition at line 42 of file MonteCarloEngine.h.

Public Member Functions

void AddChains ()
 A method to add the observable values and weights to the chain information. More...
 
void CheckHistogram (TH1 &hist, const std::string name)
 This member checks if there is overflow of the 1D histogram. More...
 
void CheckHistogram (TH2 &hist, const std::string name)
 This member checks if there is overflow of the 2D histogram. More...
 
double computeNormalizationLME ()
 A method to calculate the normalization based on the Laplace-Metropolis Estimator. More...
 
std::vector< double > computeNormalizationMC (int NIterationNormalizationMC)
 A method to calculate the normalization based on the Monte Carlo Simulation. More...
 
std::string computeStatistics ()
 A get method to compute the mean and rms of the computed observables. More...
 
void CreateHistogramMaps ()
 Creation of the histogram maps for Observables, Observable2D and Correlated Gaussian Observable. More...
 
void DefineParameters ()
 A member to classify the prior of the model parameters varied in the Monte Carlo. More...
 
double FirstDerivative (BCParameter par, std::vector< double > point)
 A method to calculate the first derivative. More...
 
double Function_h (std::vector< double > point)
 A method to calculate the LogLikelihood + LogAprioriProbability. More...
 
int getchainedObsSize () const
 A get method to access the number of Observable chains requested. More...
 
std::map< std::string, BCH1D > getHistograms1D () const
 A get method to access the map of 1D histograms. More...
 
std::map< std::string, BCH2D > getHistograms2D () const
 A get method to access the map of 2D histograms. More...
 
std::string getHistoLog () const
 A get method to access the stream that stores the log messages coming from histogram printing and checking. More...
 
int getMPIWorldSize ()
 A get method to get the number of MPI world size. More...
 
int getNumOfDiscardedEvents () const
 A get method to access the number of events discarded due to failure to update model. These events are not used for the MCMC run. More...
 
int getNumOfUsedEvents () const
 A get method to access the number of events used in the MCMC run. More...
 
bool getWriteLogLikelihoodChain ()
 A get method to get the value of the bool the writing of Loglikelihood in the ROOT tree. More...
 
bool getWriteParametersChain ()
 A get method to get the value of the bool the writing of parameters in the ROOT tree. More...
 
void InChainFillObservablesTree ()
 
void InChainFillParametersTree ()
 
void Initialize (StandardModel *Mod_i)
 Initialization of the Monte Carlo Engine. More...
 
double LogLikelihood (const std::vector< double > &parameters)
 This member calculates the loglikelihood for the observables included in the MCMC run. More...
 
void MCMCUserIterationInterface ()
 Overloaded from BCEngineMCMC in BAT More...
 
 MonteCarloEngine (const std::vector< ModelParameter > &ModPars_i, boost::ptr_vector< Observable > &Obs_i, std::vector< Observable2D > &Obs2D_i, std::vector< CorrelatedGaussianObservables > &CGO_i, std::vector< CorrelatedGaussianParameters > &CGP_i)
 Constructor. More...
 
void Print1D (BCH1D hist1D, const char *filename, int ww=0, int wh=0)
 
void Print2D (BCH2D hist2D, const char *filename, int ww=0, int wh=0)
 
void PrintCorrelationMatrixToLaTeX (const std::string filename)
 This member generates the correlation matrix using BCH2D from the BAT libraries. More...
 
void PrintHistogram (std::string &OutFile, const std::string OutputDir)
 Member used for printing histograms for observables, observable2D, correlated Gaussian observables and model parameters vs. observables. More...
 
void PrintHistogram (std::string &OutFile, Observable &it, const std::string OutputDir)
 Overloaded from PrintHistogram(BCModelOutput&, const std::string) to print histogram for observables. More...
 
double SecondDerivative (BCParameter par1, BCParameter par2, std::vector< double > point)
 A method to calculate the second derivative. More...
 
void setAlpha2D (double alpha)
 A set method to toggle the printing of legends in 1D and 2D histograms. More...
 
void setDParsFromParameters (const std::vector< double > &parameters, std::map< std::string, double > &DPars_i)
 A method to rotate the diagonalized parameters to the original basis for correlated parameters. More...
 
void setHistogram2DType (int type)
 A set method to set the band fill type for 2D histograms. More...
 
void setHistogramBufferSize (unsigned int histogramBufferSize_i)
 A set method to set the size of the buffer used by the histograms. More...
 
void setNBins1D (unsigned int nbins)
 A set method to set the number of bins for a 1D histograms. More...
 
void setNBins2D (unsigned int nbins)
 A set method to set the number of bins for a 2D histograms. More...
 
void setNChains (unsigned int i)
 A set method to fix the number of chains. More...
 
void setNoLegend (bool legend)
 A set method to toggle the printing of legends in 1D and 2D histograms. More...
 
void setPrintLoglikelihoodPlots (bool plot)
 A set method to toggle the printing of Parameters vs. Loglikelihood. More...
 
void setPrintLogo (bool print)
 A set method to toggle the printing of logo on the histogram pdf. More...
 
void setSignificants (unsigned int significants_i)
 A set method to set the number of significant digits in the output. More...
 
void setSmooth (int int_N)
 A set method to set the number of smoothing passes for ROOT in 1D histograms. More...
 
void setWriteLogLikelihoodChain (bool LL)
 A set method to toggle the writing of Loglikelihood in the ROOT tree. More...
 
void setWriteParametersChain (bool LL)
 A set method to toggle the writing of parameters in the ROOT tree. More...
 
std::string writePreRunData ()
 A method to write in a text file the best fit parameters and the prerun scale factors. More...
 
 ~MonteCarloEngine ()
 The default destructor. Some pointers defined in this class are explicitly freed. More...
 

Private Member Functions

int getPrecision (double value, double rms)
 A get method to calculate the precision of the numbers in the Statistics file based on the rms. More...
 

Private Attributes

double alpha2D
 A number between 0. and 1. that sets the opacity level of 2D Histograms, 1. being fully opaque. More...
 
std::vector< CorrelatedGaussianObservables > & CGO
 A vector of correlated Gaussian observables. More...
 
const std::vector< CorrelatedGaussianParameters > & CGP
 A vector of correlated Gaussian parameters. More...
 
unsigned int cindex
 An index to distinguish between succesive canvases used to draw histograms. More...
 
std::map< std::string, TPrincipal * > CorrelationMap
 A map between the name of a theory observable and its maximum value. More...
 
std::map< std::string, double > DPars
 A map between parameter names and their values. More...
 
std::vector< std::map< std::string, double > > DPars_allChains
 A vector of maps between parameter names and their values for all chains. More...
 
int gIdx
 
TColor * HEPfit_green
 < The number of significant digits in the Statistics File. More...
 
TColor * HEPfit_red
 < The colour green for HEPfit. More...
 
std::map< std::string, BCH1D > Histo1D
 A map between pointers to objects of type BCH1D (BAT) and their given names. More...
 
std::map< std::string, BCH2D > Histo2D
 A map between pointers to objects of type BCH2D (BAT) and their given names. More...
 
int histogram2Dtype
 Type of 2D Histogram 1001 -> box pixel, 101 -> filled, 1 -> contour. More...
 
unsigned int histogramBufferSize
 < The colour red for HEPfit. More...
 
std::ostringstream HistoLog
 A stream to store the output messages from printing and checking histograms. More...
 
double hMCMCLogLikelihood
 A variable containing the LogLikelihood values to be put into the ROOT tree. More...
 
double hMCMCLogPriorProbability
 A variable containing the LogPriorProbability values to be put into the ROOT tree. More...
 
double hMCMCLogProbability
 A variable containing the LogProbability values to be put into the ROOT tree. More...
 
std::vector< std::vector< double > > hMCMCObservables
 A vector of vectors containing the observables values of all the chains to be put into the ROOT tree. More...
 
std::vector< std::vector< double > > hMCMCObservables_weight
 A vector of vectors containing the observables weight of all the chains to be put into the ROOT tree. More...
 
TTree * hMCMCObservableTree
 A ROOT tree that contains the observables values and weight when the chains are written. More...
 
std::vector< std::vector< double > > hMCMCParameters
 A vector of vectors containing the parameter values of all the chains to be put into the ROOT tree. More...
 
TTree * hMCMCParameterTree
 A ROOT tree that contains the parameter values when the chains are written. More...
 
std::vector< double > hMCMCTree_Observables
 A vector containing the observables values to be put into the ROOT tree. More...
 
std::vector< double > hMCMCTree_Observables_weight
 A vector containing the observables weight to be put into the ROOT tree. More...
 
std::vector< double > hMCMCTree_Parameters
 A vector containing the parameter values to be put into the ROOT tree. More...
 
unsigned int kchainedObs
 The number of observables for which the chains should be written. More...
 
unsigned int kmax
 The number of observables. More...
 
unsigned int kwmax
 The number of observables whose weights are used for the MCMC. More...
 
double LogLikelihood_max
 < vector to store the value of the parameters at maximum LogLikelihood; More...
 
StandardModel * Mod
 A pointer to an abject of type Model. More...
 
const std::vector< ModelParameter > & ModPars
 A vector of model parameters. More...
 
unsigned int nBins1D
 The number of bins in a 1D histogram. More...
 
unsigned int nBins2D
 The number of bins in a 2D histogram. More...
 
bool noLegend
 A flag to toggle the histogram legends. More...
 
int nSmooth
 The number of times a 1D histogram should be smoothed. More...
 
int NumOfDiscardedEvents
 The number of events for which the update of the model fails and these events are not used for the MCMC run. More...
 
int NumOfUsedEvents
 The number of events for which the model is successfully updated and hence used for the MCMC run. More...
 
std::vector< Observable2D > & Obs2D_ALL
 A vector of all pairs of observable for Observable2D. More...
 
boost::ptr_vector< Observable > & Obs_ALL
 A vector of all observables. More...
 
double * obval
 A pointer to the vector of observable values. More...
 
double * obweight
 A pointer to the vector of observable weights. More...
 
std::ofstream ofi
 
std::vector< double > par_at_LL_max
 < The size of the buffer used for the histograms. More...
 
bool PrintLoglikelihoodPlots
 A flag to toggle the printing of Parameters vs. Loglikelihood. More...
 
bool printLogo
 A flag that is set to true for printing the logo on the histogram pdf. More...
 
int rank
 Rank of the process for a MPI run. Value is 0 for a serial run. More...
 
int rIdx
 
unsigned int significants
 
std::map< std::string, double > thMax
 A map between the name of a theory observable and its maximum value. More...
 
std::map< std::string, double > thMin
 A map between the name of a theory observable and its minimum value. More...
 
std::vector< std::string > unknownParameters
 A vector to contain the unkenown parameters passed in the configuration file. More...
 
bool WriteLogLikelihoodChain
 A flag to toggle the writing of Loglikelihood chains in the ROOT tree. More...
 
bool WriteParametersChain
 A flag to toggle the writing of parameters chains in the ROOT tree. More...
 

Constructor & Destructor Documentation

◆ MonteCarloEngine()

MonteCarloEngine::MonteCarloEngine ( const std::vector< ModelParameter > &  ModPars_i,
boost::ptr_vector< Observable > &  Obs_i,
std::vector< Observable2D > &  Obs2D_i,
std::vector< CorrelatedGaussianObservables > &  CGO_i,
std::vector< CorrelatedGaussianParameters > &  CGP_i 
)

Constructor.

Parameters
[in]ModPars_ithe vector of model parameters as defined in SomeModel.conf
[in]Obs_ithe vector of observables as defined in SomeModel.conf
[in]Obs2D_ithe vector of observables2D as defined in SomeModel.conf
[in]CGO_ithe vector of correlated Gaussian observables as defined in SomeModel.conf

Definition at line 28 of file MonteCarloEngine.cpp.

34 : BCModel(""), ModPars(ModPars_i), CGP(CGP_i), Obs_ALL(Obs_i), Obs2D_ALL(Obs2D_i),
35  CGO(CGO_i), NumOfUsedEvents(0), NumOfDiscardedEvents(0) {
36  obval = NULL;
37  obweight = NULL;
38  Mod = NULL;
39  cindex = 0;
40  printLogo = false;
41  nSmooth = 0;
42  histogram2Dtype = 1001;
43  noLegend = true;
44  PrintLoglikelihoodPlots = false;
45  WriteLogLikelihoodChain = false;
46  WriteParametersChain = false;
47  alpha2D = 1.;
48  kchainedObs = 0;
49  nBins1D = NBINS1D;
50  nBins2D = NBINS2D;
51  significants = 0;
52  histogramBufferSize = 0;
53  LogLikelihood_max = std::numeric_limits<double>::lowest();
54 #ifdef _MPI
55  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
56 #else
57  rank = 0;
58 #endif
59  if (rank == 0) {
60  TH1::StatOverflows(kTRUE);
61  TH1::SetDefaultBufferSize(100000);
62 
63 #if ROOT_VERSION_CODE > ROOT_VERSION(6,0,0)
64  gIdx = TColor::GetFreeColorIndex();
65  rIdx = TColor::GetFreeColorIndex() + 1;
66 #else
67  gIdx = 1000;
68  rIdx = 1001;
69 #endif
70 
71  HEPfit_green = new TColor(gIdx, 0.0, 0.56, 0.57, "HEPfit_green");
72  HEPfit_red = new TColor(rIdx, 0.57, 0.01, 0.00, "HEPfit_red");
73  }
74 };

◆ ~MonteCarloEngine()

MonteCarloEngine::~MonteCarloEngine ( )

The default destructor. Some pointers defined in this class are explicitly freed.

Definition at line 186 of file MonteCarloEngine.cpp.

187 {
188  if (rank == 0) { // These are created only by the master
189  delete [] obval;
190  delete [] obweight;
191  delete HEPfit_red;
192  delete HEPfit_green;
193  HEPfit_red = NULL;
194  HEPfit_green = NULL;
195  if (CorrelationMap.size() > 0) {
196  for (std::map<std::string, TPrincipal *>::iterator it = CorrelationMap.begin(); it != CorrelationMap.end(); it++) {
197  delete it->second;
198  it->second = NULL;
199  }
200  }
201  }
202 };

Member Function Documentation

◆ AddChains()

void MonteCarloEngine::AddChains ( )

A method to add the observable values and weights to the chain information.

If the WriteChain flag is set to true in the MonteCarlo.conf file then chain information is written. This method adds the observable values and weights (for the observables included in the MCMC) to the chain information.

Definition at line 884 of file MonteCarloEngine.cpp.

884  {
885  if (fMCMCFlagWriteChainToFile) InitializeMarkovChainTree();
886  TDirectory* dir = gDirectory;
887  GetOutputFile()->cd();
888 
889  hMCMCObservableTree = new TTree(TString::Format("%s_Observables", GetSafeName().data()), TString::Format("%s_Observables", GetSafeName().data()));
890  hMCMCObservableTree->Branch("Chain", &fMCMCTree_Chain, "chain/i");
891  hMCMCObservableTree->Branch("Iteration", &fMCMCCurrentIteration, "iteration/i");
892  if (WriteLogLikelihoodChain) {
893  hMCMCObservableTree->Branch("LogLikelihood", &hMCMCLogLikelihood, "loglikelihood/D");
894  hMCMCObservableTree->Branch("LogProbability", &hMCMCLogProbability, "logprobability/D");
895  hMCMCObservableTree->Branch("LogPriorProbability", &hMCMCLogPriorProbability, "logpriorprobability/D");
896  }
897  if (fMCMCFlagWriteChainToFile) {
898  hMCMCObservables.assign(fMCMCNChains, std::vector<double>(Obs_ALL.size(), 0.));
899  hMCMCTree_Observables.assign(Obs_ALL.size(), 0.);
900  if (kwmax > 0) {
901  hMCMCObservables_weight.assign(fMCMCNChains, std::vector<double>(kwmax, 0.));
902  hMCMCTree_Observables_weight.assign(kwmax, 0.);
903  }
904  int k = 0, kweight = 0;
905  for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin(); it < Obs_ALL.end(); it++) {
906  hMCMCObservableTree->Branch(it->getName().data(), &hMCMCTree_Observables[k], (it->getName() + "/D").data());
907  hMCMCObservableTree->SetAlias(TString::Format("HEPfit_Observables%i", k), it->getName().data());
908  k++;
909  if (!it->isTMCMC() && it->getDistr().compare("weight") == 0) {
910  hMCMCObservableTree->Branch((it->getName() + "_weight").data(), &hMCMCTree_Observables_weight[kweight], (it->getName() + "_weight/D").data());
911  hMCMCObservableTree->SetAlias(TString::Format("HEPfit_Observables_weight%i", kweight), (it->getName() + "_weight").data());
912  kweight++;
913  }
914  }
915  } else if (getchainedObsSize() > 0) {
916  hMCMCObservables.assign(fMCMCNChains, std::vector<double>(getchainedObsSize(), 0.));
917  hMCMCTree_Observables.assign(getchainedObsSize(), 0.);
918  int k = 0;
919  for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin(); it < Obs_ALL.end(); it++) {
920  if (it->isWriteChain()) {
921  hMCMCObservableTree->Branch(it->getName().data(), &hMCMCTree_Observables[k], (it->getName() + "/D").data());
922  hMCMCObservableTree->SetAlias(TString::Format("HEPfit_Observables%i", k), it->getName().data());
923  k++;
924  }
925  }
926  }
927  hMCMCObservableTree->SetAutoSave(10 * fMCMCNIterationsPreRunCheck);
928  hMCMCObservableTree->AutoSave("SelfSave");
929 
930  hMCMCParameterTree = new TTree(TString::Format("%s_Parameters", GetSafeName().data()), TString::Format("%s_Parameters", GetSafeName().data()));
931  hMCMCParameterTree->Branch("Chain", &fMCMCTree_Chain, "chain/i");
932  hMCMCParameterTree->Branch("Iteration", &fMCMCCurrentIteration, "iteration/i");
933  if (WriteParametersChain) {
934  hMCMCParameters.assign(fMCMCNChains, std::vector<double>(DPars.size(), 0.));
935  hMCMCTree_Parameters.assign(DPars.size(), 0.);
936  int k = 0;
937  for (std::map<std::string, double>::iterator it = DPars.begin(); it != DPars.end(); it++) {
938  hMCMCParameterTree->Branch(it->first.data(), &hMCMCTree_Parameters[k], (it->first + "/D").data());
939  hMCMCParameterTree->SetAlias(TString::Format("HEPfit_Parameters%i", k), it->first.data());
940  k++;
941  }
942  }
943  hMCMCParameterTree->SetAutoSave(10 * fMCMCNIterationsPreRunCheck);
944  hMCMCParameterTree->AutoSave("SelfSave");
945 
946  gDirectory = dir;
947 }

◆ CheckHistogram() [1/2]

void MonteCarloEngine::CheckHistogram ( TH1 &  hist,
const std::string  name 
)

This member checks if there is overflow of the 1D histogram.

Parameters
[in]hista reference to an object of type TH1D as defined in the ROOT libraries
[in]namethe name for the histogram

Definition at line 594 of file MonteCarloEngine.cpp.

594  {
595  double UnderFlowContent = hist.GetBinContent(0);
596  double OverFlowContent = hist.GetBinContent(nBins1D + 1);
597  double Integral = hist.Integral();
598  double TotalContent = 0.0;
599  for (unsigned int n = 0; n <= nBins1D + 1; n++)
600  TotalContent += hist.GetBinContent(n);
601  HistoLog << name << ": "
602  << Integral / TotalContent * 100. << "% within the range, "
603  << UnderFlowContent / TotalContent * 100. << "% underflow, "
604  << OverFlowContent / TotalContent * 100. << "% overflow"
605  << std::endl;
606 }

◆ CheckHistogram() [2/2]

void MonteCarloEngine::CheckHistogram ( TH2 &  hist,
const std::string  name 
)

This member checks if there is overflow of the 2D histogram.

Parameters
[in]hista reference to an object of type TH2D as defined in the ROOT libraries
[in]namethe name for the histogram

Definition at line 608 of file MonteCarloEngine.cpp.

608  {
609  double Integral = hist.Integral();
610  double TotalContent = 0.0;
611  for (unsigned int m = 0; m <= nBins2D + 1; m++)
612  for (unsigned int n = 0; n <= nBins2D + 1; n++)
613  TotalContent += hist.GetBinContent(m, n);
614  HistoLog << name << ": "
615  << Integral / TotalContent * 100. << "% within the ranges"
616  << std::endl;
617 }

◆ computeNormalizationLME()

double MonteCarloEngine::computeNormalizationLME ( )

A method to calculate the normalization based on the Laplace-Metropolis Estimator.

Returns
the normalization

Definition at line 1335 of file MonteCarloEngine.cpp.

1335  {
1336 /* PENDING REVIEW FOR USE WITH BAT v1.0. */
1337  unsigned int Npars = GetNParameters();
1338  std::vector<double> mode(GetBestFitParameters());
1339  if (mode.size() == 0) {
1340  throw std::runtime_error("\n ERROR: Global Mode could not be determined possibly because of infinite loglikelihood. Normalization computation cannot be completed.\n");
1341  }
1342  gslpp::matrix<double> Hessian(Npars, Npars, 0.);
1343 
1344  for (unsigned int i = 0; i < Npars; i++)
1345  for (unsigned int j = 0; j < Npars; j++) {
1346  // calculate Hessian matrix element
1347  Hessian.assign(i, j, -SecondDerivative(GetParameter(i), GetParameter(j), mode));
1348  }
1349  double det_Hessian = Hessian.determinant();
1350 
1351  return exp(Npars / 2. * log(2. * M_PI) + 0.5 * log(1. / det_Hessian) + LogLikelihood(mode) + LogAPrioriProbability(mode));
1352 }

◆ computeNormalizationMC()

std::vector< double > MonteCarloEngine::computeNormalizationMC ( int  NIterationNormalizationMC)

A method to calculate the normalization based on the Monte Carlo Simulation.

Returns
the normalization

Definition at line 1319 of file MonteCarloEngine.cpp.

1319  {
1320  // Number of MC iterations
1321  SetNIterationsMin(NIterationNormalizationMC);
1322  SetIntegrationMethod(BCIntegrate::kIntMonteCarlo);
1323  Integrate();
1324  std::vector<double> norm;
1325  norm.clear();
1326  norm.push_back(GetIntegral());
1327  norm.push_back(GetError());
1328  if (norm[0] < 0.) {
1329  throw std::runtime_error("\n ERROR: Normalization computation cannot be completed since integral is negative.\n");
1330  }
1331 
1332  return norm;
1333 }

◆ computeStatistics()

std::string MonteCarloEngine::computeStatistics ( )

A get method to compute the mean and rms of the computed observables.

Returns
a string containing the statistics

Definition at line 1012 of file MonteCarloEngine.cpp.

1012  {
1013 
1014  std::vector<double> mode(GetBestFitParameters());
1015  if (mode.size() == 0) throw std::runtime_error("\n ERROR: Global Mode could not be determined possibly because of infinite loglikelihood. Observables statistics cannot be generated.\n");
1016  std::streamsize ss_prec = std::cout.precision();
1017  unsigned int rmsPrecision = 2;
1018  if (significants > 0) rmsPrecision = significants;
1019  std::ostringstream StatsLog;
1020  int i = 0;
1021  StatsLog << "Statistics file for Observables, Binned Observables and Correlated Gaussian Observables.\n" << std::endl;
1022  if (Obs_ALL.size() > 0) StatsLog << "Observables:\n" << std::endl;
1023  for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin(); it < Obs_ALL.end(); it++) {
1024  StatsLog.precision(ss_prec); /* resets precision*/
1025  if (it->getObsType().compare("BinnedObservable") == 0) {
1026  StatsLog << " (" << ++i << ") Binned Observable \"";
1027  StatsLog << it->getName() << "[" << it->getTho()->getBinMin() << ", " << it->getTho()->getBinMax() << "]" << "\":";
1028  } else if (it->getObsType().compare("FunctionObservable") == 0) {
1029  StatsLog << " (" << ++i << ") Function Observable \"";
1030  StatsLog << it->getName() << "[" << it->getTho()->getBinMin() << "]" << "\":";
1031  } else if (it->getObsType().compare("HiggsObservable") == 0) {
1032  StatsLog << " (" << ++i << ") Higgs Observable \"";
1033  StatsLog << it->getName() << "\":";
1034  } else {
1035  StatsLog << " (" << ++i << ") " << it->getObsType() << " \"";
1036  StatsLog << it->getName() << "\":";
1037  }
1038  StatsLog << std::endl;
1039 
1040  BCH1D bch1d = Histo1D[it->getName()];
1041 
1042  if (bch1d.GetHistogram()->Integral() > 0.0) {
1043  double rms = bch1d.GetHistogram()->GetRMS();
1044  StatsLog << " Mean +- sqrt(V): " << std::setprecision(getPrecision(bch1d.GetHistogram()->GetMean(),rms))
1045  << bch1d.GetHistogram()->GetMean() << " +- " << std::setprecision(rmsPrecision)
1046  << rms << std::endl
1047 
1048  << " (Marginalized) mode: " << std::setprecision(getPrecision(bch1d.GetLocalMode(0),rms)) << bch1d.GetLocalMode(0) << std::endl;
1049  std::vector<double> intervals;
1050  intervals.push_back(0.682689492137);
1051  intervals.push_back(0.954499736104);
1052  intervals.push_back(0.997300203937);
1053  std::vector<BCH1D::BCH1DSmallestInterval> v = bch1d.GetSmallestIntervals(intervals);
1054  for (unsigned int i = 0; i < v.size(); i++) {
1055  StatsLog << " Smallest interval(s) containing at least " << std::setprecision(ss_prec) << v[i].total_mass * 100 << "% and local mode(s):"
1056  << std::endl;
1057  for (unsigned j = 0; j < v[i].intervals.size(); j++) {
1058  double interval_xmin = v[i].intervals[j].xmin;
1059  double interval_xmax = v[i].intervals[j].xmax;
1060  double interval_mode = v[i].intervals[j].mode;
1061  double interval_heignt = v[i].intervals[j].relative_height;
1062  double interval_relative_mass = v[i].intervals[j].relative_mass;
1063  StatsLog << " (" << std::setprecision(getPrecision(interval_xmin, rms)) << interval_xmin << ", " << std::setprecision(getPrecision(interval_xmax, rms)) << interval_xmax
1064  << ") (local mode at " << std::setprecision(getPrecision(interval_mode, rms)) << interval_mode << " with rel. height "
1065  << std::setprecision(getPrecision(interval_heignt, ss_prec)) << interval_heignt << "; rel. area " << std::setprecision(getPrecision(interval_relative_mass, ss_prec)) << interval_relative_mass << ")"
1066  << std::endl;
1067  StatsLog << std::endl;
1068  }
1069  }
1070  } else {
1071  StatsLog << "\nWARNING: The histogram of " << it->getName() << " is empty! Statistics cannot be generated\n" << std::endl;
1072  }
1073  }
1074 
1075  if (CGO.size() > 0) StatsLog << "\nCorrelated (Gaussian) Observables:\n" << std::endl;
1076  for (std::vector<CorrelatedGaussianObservables>::iterator it1 = CGO.begin(); it1 < CGO.end(); it1++) {
1077  StatsLog << "\n" << it1->getName() << ":\n" << std::endl;
1078  i = 0;
1079  std::vector<Observable> CGObs(it1->getObs());
1080  for (std::vector<Observable>::iterator it2 = CGObs.begin(); it2 < CGObs.end(); it2++) {
1081  StatsLog.precision(ss_prec); /* resets precision*/
1082  if (it2->getObsType().compare("BinnedObservable") == 0) {
1083  StatsLog << " (" << ++i << ") Binned Observable \"";
1084  StatsLog << it2->getName() << "[" << it2->getTho()->getBinMin() << ", " << it2->getTho()->getBinMax() << "]" << "\":";
1085  } else if (it2->getObsType().compare("FunctionObservable") == 0) {
1086  StatsLog << " (" << ++i << ") Function Observable \"";
1087  StatsLog << it2->getName() << "[" << it2->getTho()->getBinMin() << "]" << "\":";
1088  } else if (it2->getObsType().compare("HiggsObservable") == 0) {
1089  StatsLog << " (" << ++i << ") Higgs Observable \"";
1090  StatsLog << it2->getName() << "\":";
1091  } else {
1092  StatsLog << " (" << ++i << ") " << it2->getObsType() << " \"";
1093  StatsLog << it2->getName() << "\":";
1094  }
1095 
1096  StatsLog << std::endl;
1097  BCH1D bch1d = Histo1D[it2->getName()];
1098  if (bch1d.GetHistogram()->Integral() > 0.0) {
1099  double rms = bch1d.GetHistogram()->GetRMS();
1100  StatsLog << " Mean +- sqrt(V): " << std::setprecision(getPrecision(bch1d.GetHistogram()->GetMean(), rms))
1101  << bch1d.GetHistogram()->GetMean() << " +- " << std::setprecision(rmsPrecision)
1102  << rms << std::endl
1103 
1104  << " (Marginalized) mode: " << std::setprecision(getPrecision(bch1d.GetLocalMode(0), rms)) << bch1d.GetLocalMode(0) << std::endl;
1105 
1106  std::vector<double> intervals;
1107  intervals.push_back(0.682689492137);
1108  intervals.push_back(0.954499736104);
1109  intervals.push_back(0.997300203937);
1110 
1111  std::vector<BCH1D::BCH1DSmallestInterval> v = bch1d.GetSmallestIntervals(intervals);
1112  for (unsigned int i = 0; i < v.size(); i++) {
1113  StatsLog << " Smallest interval(s) containing at least " << std::setprecision(ss_prec) << v[i].total_mass * 100 << "% and local mode(s):" << std::endl;
1114  for (unsigned j = 0; j < v[i].intervals.size(); j++) {
1115  double interval_xmin = v[i].intervals[j].xmin;
1116  double interval_xmax = v[i].intervals[j].xmax;
1117  double interval_mode = v[i].intervals[j].mode;
1118  double interval_heignt = v[i].intervals[j].relative_height;
1119  double interval_relative_mass = v[i].intervals[j].relative_mass;
1120  StatsLog << " (" << std::setprecision(getPrecision(interval_xmin, rms)) << interval_xmin << ", " << std::setprecision(getPrecision(interval_xmax, rms)) << interval_xmax
1121  << ") (local mode at " << std::setprecision(getPrecision(interval_mode, rms)) << interval_mode << " with rel. height "
1122  << std::setprecision(getPrecision(interval_heignt, ss_prec)) << interval_heignt << "; rel. area " << std::setprecision(getPrecision(interval_relative_mass, ss_prec)) << interval_relative_mass << ")"
1123  << std::endl;
1124  StatsLog << std::endl;
1125  }
1126  }
1127  } else {
1128  StatsLog << "\nWARNING: The histogram of " << it2->getName() << " is empty! Statistics cannot be generated\n" << std::endl;
1129  }
1130  }
1131  if (it1->isPrediction()) {
1132  int size = it1->getObs().size();
1133  CorrelationMap[it1->getName()]->MakePrincipals();
1134  //CorrelationMap[it1->getName()]->Print();
1135  TMatrixD * corr = const_cast<TMatrixD*>(CorrelationMap[it1->getName()]->GetCovarianceMatrix()); // This returns the normalized correlation matrix.
1136  TVectorD * sigma = const_cast<TVectorD*>(CorrelationMap[it1->getName()]->GetSigmas()); // This returns the vector of standard deviations.
1137  *corr *= (double)size; // Get rid of the normalization which is just the size of the matrix.
1138  gslpp::matrix<double> inverseCovariance(size, size);
1139  for (int i = 0; i < size; i++) {
1140  for (int j = 0; j <= i; j++) {
1141  inverseCovariance(i, j) = (*corr)(i, j) * (*sigma)(i) * (*sigma)(j);
1142  inverseCovariance(j, i) = inverseCovariance(i, j);
1143  }
1144  }
1145  bool SingularCovariance = inverseCovariance.isSingular();
1146  if (!SingularCovariance) inverseCovariance = inverseCovariance.inverse(); // This is just produces inverse covariance, the name is misleading.
1147  StatsLog << "\nThe correlation matrix for " << it1->getName() << " is given by the " << size << "x"<< size << " matrix:\n" << std::endl;
1148 
1149  for (int i = 0; i < size + 1; i++) {
1150  if (i == 0) StatsLog << std::setw(4) << "" << " | ";
1151  else StatsLog << std::setw(6) << i << std::setw(6) << " |";
1152  }
1153  StatsLog << std::endl;
1154  for (int i = 0; i < size + 1; i++) {
1155  if (i == 0) StatsLog << std::setw(8) << "--------";
1156  else StatsLog << std::setw(12) << "------------";
1157  }
1158  StatsLog << std::endl;
1159  for (int i = 0; i < size; i++) {
1160  for (int j = 0; j < size + 1; j++) {
1161  if (j == 0) StatsLog << std::setw(4) << i+1 << " |";
1162  else StatsLog << std::setprecision(5) << std::setw(12) << (*corr)(i, j - 1);
1163  }
1164  StatsLog << std::endl;
1165  }
1166  StatsLog << std::endl;
1167  if (!SingularCovariance) {
1168  StatsLog << " The inverse of the square root of the diagonal elements of the inverse covariance matrix are:\n" << std::endl;
1169  for (int i = 0; i < size + 1; i++) {
1170  if (i == 0) StatsLog << std::setw(4) << "sigma" << "|";
1171  else StatsLog << std::setprecision(5) << std::setw(12) << 1. / sqrt(inverseCovariance(i - 1, i - 1));
1172  }
1173  } else StatsLog << " The covariance matrix cannot be inverted.\n" << std::endl;
1174  StatsLog << std::endl;
1175  }
1176  }
1177 
1178  setDParsFromParameters(mode,DPars);
1179  Mod->Update(DPars);
1180 
1181  // values for controlling format
1182  int par_width = 0;
1183  int obs_width = 0;
1184  int value_width = 15;
1185  for (std::map<std::string,double>::iterator it = DPars.begin(); it != DPars.end(); it++) par_width = std::max(par_width, (int)(it->first).size());
1186  for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin(); it < Obs_ALL.end(); it++) obs_width = std::max(obs_width, (int)(it->getName()).size());
1187  par_width = par_width + 1;
1188  obs_width = obs_width + 1;
1189  const std::string sep = " |";
1190  const std::string par_line = sep + std::string(par_width + value_width + sep.size() * 2 - 1, '-') + '|';
1191  const std::string obs_line = sep + std::string(obs_width + value_width + sep.size() * 2 - 1, '-') + '|';
1192  StatsLog << std::setprecision(5);
1193 
1194  StatsLog << "\n*** Statistical details using global mode ***\n" << std::endl;
1195 
1196  StatsLog << "\nValue of the parameters at the global mode:" << std::endl;
1197  StatsLog << std::endl;
1198  StatsLog << par_line << '\n' << sep
1199  << std::left << std::setw(par_width) << "parameter" << sep << std::right << std::setw(value_width) << "value at mode" << sep << '\n' << par_line << '\n';
1200 
1201  for (std::map<std::string,double>::iterator it = DPars.begin(); it != DPars.end(); it++)
1202  StatsLog << sep << std::left << std::setw(par_width) << it->first << sep << std::right << std::setw(value_width) << it->second << sep << '\n';
1203 
1204  StatsLog << par_line << '\n';
1205  StatsLog << std::endl;
1206 
1207  StatsLog << "\nValue of the observables at the global mode:" << std::endl;
1208  StatsLog << std::endl;
1209  StatsLog << obs_line << '\n' << sep
1210  << std::left << std::setw(obs_width) << "observable" << sep << std::right << std::setw(value_width) << "value at mode" << sep << '\n' << obs_line << '\n';
1211 
1212  for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin(); it < Obs_ALL.end(); it++)
1213  StatsLog << sep << std::left << std::setw(obs_width) << it->getName() << sep << std::right << std::setw(value_width) << it->computeTheoryValue() << sep << '\n';
1214 
1215  StatsLog << obs_line << '\n';
1216  StatsLog << std::endl;
1217 
1218  StatsLog << "LogProbability at mode: " << LogLikelihood(mode) + LogAPrioriProbability(mode) << std::endl;
1219  StatsLog << "LogLikelihood at mode: " << LogLikelihood(mode) << std::endl;
1220  StatsLog << "LogAPrioriProbability at mode: " << LogAPrioriProbability(mode) << "\n\n" << std::endl;
1221 
1222  double llika = Histo1D["LogLikelihood"].GetHistogram()->GetMean();
1223  StatsLog << "LogLikelihood mean value: " << llika << std::endl;
1224  double llikv = Histo1D["LogLikelihood"].GetHistogram()->GetRMS();
1225  llikv *= llikv;
1226  StatsLog << "LogLikelihood variance: " << llikv << std::endl;
1227  double dbar = -2.*llika; //Wikipedia notation...
1228  double pd = 2.*llikv; //Wikipedia notation...
1229  StatsLog << "IC value: " << dbar + 2.*pd << std::endl;
1230  StatsLog << "DIC value: " << dbar + pd << std::endl;
1231  StatsLog << std::endl;
1232  StatsLog << std::endl;
1233 
1234 
1235  //For testing purposes:
1236  const BCEngineMCMC::Statistics& st = GetStatistics();
1237  //get mean value of parameters from BAT
1238  std::vector<double> parmeans = st.mean;
1239 
1240  setDParsFromParameters(parmeans,DPars);
1241  Mod->Update(DPars);
1242 
1243  StatsLog << "*** Statistical details using mean values of parameters ***\n" << std::endl;
1244 
1245  StatsLog << "\nMean value of the parameters:" << std::endl;
1246  StatsLog << std::endl;
1247  StatsLog << par_line << '\n' << sep
1248  << std::left << std::setw(par_width) << "parameter" << sep << std::right << std::setw(value_width) << "mean value" << sep << '\n' << par_line << '\n';
1249 
1250  for (std::map<std::string,double>::iterator it = DPars.begin(); it != DPars.end(); it++)
1251  StatsLog << sep << std::left << std::setw(par_width) << it->first << sep << std::right << std::setw(value_width) << it->second << sep << '\n';
1252 
1253  StatsLog << par_line << '\n';
1254  StatsLog << std::endl;
1255 
1256  StatsLog << "Mean of LogProbability: " << st.probability_mean << std::endl;
1257  StatsLog << "Variance of LogProbability: " << st.probability_variance << std::endl;
1258  StatsLog << "LogProbability at mode: " << st.probability_at_mode << std::endl;
1259  StatsLog << std::endl;
1260 
1261  double llonmean = LogLikelihood(parmeans);
1262  StatsLog << "LogLikelihood on mean value of parameters: " << llonmean << std::endl;
1263  StatsLog << "pD computed using variance: " << pd << std::endl;
1264  pd = 2.*llonmean-2.*llika;
1265  StatsLog << "pD computed using 2LL(thetabar) - 2LLbar: " << pd << std::endl;
1266  StatsLog << "IC value computed from BAT with alternate pD definition: " << dbar + 2.*pd << std::endl;
1267  StatsLog << "DIC value computed from BAT with alternate pD definition: " << dbar + pd << std::endl;
1268  StatsLog << std::endl;
1269  StatsLog << std::endl;
1270 
1271 
1272  setDParsFromParameters(par_at_LL_max,DPars);
1273  Mod->Update(DPars);
1274 
1275  StatsLog << "*** Statistical details using parameter values at maximum LogLikelihood ***\n" << std::endl;
1276 
1277  StatsLog << "\nValue of the parameters at maximum LogLikelihood:" << std::endl;
1278  StatsLog << std::endl;
1279  StatsLog << par_line << '\n' << sep
1280  << std::left << std::setw(par_width) << "parameter" << sep << std::right << std::setw(value_width) << "value at max." << sep << '\n' << par_line << '\n';
1281 
1282  for (std::map<std::string,double>::iterator it = DPars.begin(); it != DPars.end(); it++)
1283  StatsLog << sep << std::left << std::setw(par_width) << it->first << sep << std::right << std::setw(value_width) << it->second << sep << '\n';
1284 
1285  StatsLog << par_line << '\n';
1286  StatsLog << std::endl;
1287 
1288  StatsLog << "\nValue of the observables at the maximum LogLikelihood:" << std::endl;
1289  StatsLog << std::endl;
1290  StatsLog << obs_line << '\n' << sep
1291  << std::left << std::setw(obs_width) << "observable" << sep << std::right << std::setw(value_width) << "value at max." << sep << '\n' << obs_line << '\n';
1292 
1293  for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin(); it < Obs_ALL.end(); it++)
1294  StatsLog << sep << std::left << std::setw(obs_width) << it->getName() << sep << std::right << std::setw(value_width) << it->computeTheoryValue() << sep << '\n';
1295 
1296  StatsLog << obs_line << '\n';
1297  StatsLog << std::endl;
1298 
1299  StatsLog << "Maximum LogLikelihood: " << LogLikelihood_max << std::endl;
1300 
1301  StatsLog << std::endl;
1302 
1303  return StatsLog.str().c_str();
1304 }

◆ CreateHistogramMaps()

void MonteCarloEngine::CreateHistogramMaps ( )

Creation of the histogram maps for Observables, Observable2D and Correlated Gaussian Observable.

Definition at line 121 of file MonteCarloEngine.cpp.

122 {
123  if (histogramBufferSize != 0) TH1::SetDefaultBufferSize(histogramBufferSize);
124 
125  TH1D * lhisto = new TH1D("LogLikelihood", "LogLikelihood", nBins1D, 1., -1.);
126  lhisto->GetXaxis()->SetTitle("LogLikelihood");
127  BCH1D bclhisto = BCH1D(lhisto);
128  Histo1D["LogLikelihood"] = bclhisto;
129 
130  for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin(); it != Obs_ALL.end(); it++) {
131  std::string HistName = it->getName();
132  if (Histo1D.find(HistName) == Histo1D.end()) {
133  TH1D * histo = new TH1D(HistName.c_str(), it->getLabel().c_str(), nBins1D, it->getMin(), it->getMax());
134  histo->GetXaxis()->SetTitle(it->getLabel().c_str());
135  BCH1D bchisto = BCH1D(histo);
136  Histo1D[HistName] = bchisto;
137  }
138  }
139  for (std::vector<Observable2D>::iterator it = Obs2D_ALL.begin(); it != Obs2D_ALL.end(); it++) {
140  std::string HistName = it->getName();
141  if (Histo2D.find(HistName) == Histo2D.end()) {
142  TH2D * histo2 = new TH2D(HistName.c_str(), (it->getLabel() + " vs. " + it->getLabel2()).c_str(), nBins2D, it->getMin(), it->getMax(), nBins2D, it->getMin2(), it->getMax2());
143  histo2->GetXaxis()->SetTitle(it->getLabel().c_str());
144  histo2->GetYaxis()->SetTitle(it->getLabel2().c_str());
145  BCH2D bchisto2 = BCH2D(histo2);
146  Histo2D[HistName] = bchisto2;
147  }
148  }
149  for (std::vector<CorrelatedGaussianObservables>::iterator it1 = CGO.begin(); it1 != CGO.end(); ++it1) {
150  std::vector<Observable> ObsV(it1->getObs());
151  for (std::vector<Observable>::iterator it = ObsV.begin(); it != ObsV.end(); ++it) {
152  std::string HistName = it->getName();
153  if (Histo1D.find(HistName) == Histo1D.end()) {
154  TH1D * histo = new TH1D(HistName.c_str(), it->getLabel().c_str(), nBins1D, it->getMin(), it->getMax());
155  histo->GetXaxis()->SetTitle(it->getLabel().c_str());
156  BCH1D bchisto = BCH1D(histo);
157  Histo1D[HistName] = bchisto;
158  }
159  }
160  }
161 
162  if (PrintLoglikelihoodPlots) {
163  for (std::vector<ModelParameter>::const_iterator it = ModPars.begin(); it != ModPars.end(); it++) {
164  if (it->IsFixed()) continue;
165  if (std::find(unknownParameters.begin(), unknownParameters.end(), it->getname()) != unknownParameters.end()) continue;
166  std::string HistName = it->getname() + "_vs_LogLikelihood";
167  if (Histo2D.find(HistName) == Histo2D.end()) {
168  TH2D * histo2 = new TH2D(HistName.c_str(), (it->getname() + " vs. LogLikelihood").c_str(), nBins2D, 1., -1., nBins2D, 1., -1.);
169  histo2->GetXaxis()->SetTitle(it->getname().c_str());
170  histo2->GetYaxis()->SetTitle("LogLikelihood");
171  BCH2D bchisto2 = BCH2D(histo2);
172  Histo2D[HistName] = bchisto2;
173  }
174  }
175  }
176 };

◆ DefineParameters()

void MonteCarloEngine::DefineParameters ( )

A member to classify the prior of the model parameters varied in the Monte Carlo.

The model parameters being varied are first sorted out checking for the existence of of Gaussian error ( \(\delta_{g}\)) or flat error ( \(\delta_{f}\)). The SetPriorGaus method and SetPriorConstant method in BCEngineMCMC (BAT) is used to set a Gaussian or flat prior for the model parameters respectively. In case a model parameter has both a Gaussian and a flat error then the combined function is built using TF1. The definition of the combined error is

\[ \delta_{combined}(x) = \frac{1}{4\delta_f}\left(\frac{Erf(x - x_{ave}+\delta_f)}{\sqrt{2}\delta_g} - \frac{Erf(x - x_{ave}-\delta_f)}{\sqrt{2}\delta_g}\right), \]

with

\[ Erf(x)=\frac{1}{\sqrt{2}\pi}\int_0^x \exp^{-\frac{t^2}{2}}dt. \]

Definition at line 206 of file MonteCarloEngine.cpp.

206  {
207  // Add parameters to your model here.
208  // You can then use them in the methods below by calling the
209  // parameters.at(i) or parameters[i], where i is the index
210  // of the parameter. The indices increase from 0 according to the
211  // order of adding the parameters.
212  if (rank == 0) std::cout << "\nParameters varied in this run:" << std::endl;
213  unsigned int k = 0;
214  for (std::vector<ModelParameter>::const_iterator it = ModPars.begin();
215  it < ModPars.end(); it++) {
216  if (std::find(unknownParameters.begin(), unknownParameters.end(), it->getname()) == unknownParameters.end()) {
217  if (it->geterrf() == 0. && it->geterrg() == 0.)
218  continue;
219 
220  AddParameter(it->getname().c_str(), it->getmin(), it->getmax());
221  if (rank == 0) std::cout << k << ": " << it->getname() << ", ";
222 
223  if (it->IsCorrelated()) {
224  for (unsigned int i = 0; i < CGP.size(); i++) {
225  if (CGP[i].getName().compare(it->getCgp_name()) == 0) {
226  std::string index = it->getname().substr(CGP[i].getName().size());
227  long int lindex = strtol(index.c_str(), NULL, 10);
228  if (lindex > 0)
229  DPars[CGP[i].getPar(lindex - 1).getname()] = 0.;
230  else {
231  std::stringstream out;
232  out << it->getname();
233  throw std::runtime_error("MonteCarloEngine::DefineParameters(): " + out.str() + "seems to be part of a CorrelatedGaussianParameters object, but I couldn't find the corresponding object");
234  }
235  }
236  }
237  } else
238  DPars[it->getname()] = 0.;
239  if (it->geterrf() == 0.) GetParameter(k).SetPrior(new BCGaussianPrior(it->getave(), it->geterrg())); //SetPriorGauss(k, it->getave(), it->geterrg());
240  else if (it->geterrg() == 0.) GetParameter(k).SetPriorConstant(); //SetPriorConstant(k);
241  else {
242  TF1 combined = TF1(it->getname().c_str(),
243  "(TMath::Erf((x-[0]+[2])/sqrt(2.)/[1])-TMath::Erf((x-[0]-[2])/sqrt(2.)/[1]))/4./[2]",
244  it->getmin(), it->getmax());
245  combined.SetParameter(0, it->getave());
246  combined.SetParameter(1, it->geterrg());
247  combined.SetParameter(2, it->geterrf());
248  GetParameter(k).SetPrior(new BCTF1Prior(combined)); //SetPrior(k, combined);
249  }
250  k++;
251  }
252 
253  }
254  if (unknownParameters.size() > 0 && rank == 0) {
255  std::cout << "\n" << std::endl;
256  for (std::vector<std::string>::iterator it = unknownParameters.begin(); it != unknownParameters.end(); it++)
257  std::cout << "WARNING: unknown parameter " << *it << " not added to MCMC" << std::endl;
258  }
259 }

◆ FirstDerivative()

double MonteCarloEngine::FirstDerivative ( BCParameter  par,
std::vector< double >  point 
)

A method to calculate the first derivative.

Returns
the first derivative

Definition at line 1389 of file MonteCarloEngine.cpp.

1389  {
1390 
1391  if (point.size() != GetNParameters()) {
1392  throw std::runtime_error("MCMCENgine::FirstDerivative : Invalid number of entries in the vector.");
1393  }
1394 
1395  // define steps
1396  const double dx1 = par.GetRangeWidth() / NSTEPS;
1397  const double dx2 = dx1 * 2.;
1398  const double dx3 = dx1 * 3.;
1399 
1400  // define points at which to evaluate
1401  std::vector<double> x1p = point;
1402  std::vector<double> x1m = point;
1403  std::vector<double> x2p = point;
1404  std::vector<double> x2m = point;
1405  std::vector<double> x3p = point;
1406  std::vector<double> x3m = point;
1407 
1408  unsigned idx = GetParameters().Index(par.GetName());
1409 
1410  x1p[idx] += dx1;
1411  x1m[idx] -= dx1;
1412  x2p[idx] += dx2;
1413  x2m[idx] -= dx2;
1414  x3p[idx] += dx3;
1415  x3m[idx] -= dx3;
1416 
1417  const double m1 = (Function_h(x1p) - Function_h(x1m)) / 2. / dx1;
1418  const double m2 = (Function_h(x2p) - Function_h(x2m)) / 4. / dx1;
1419  const double m3 = (Function_h(x3p) - Function_h(x3m)) / 6. / dx1;
1420 
1421  return 3. / 2. * m1 - 3. / 5. * m2 + 1. / 10. * m3;
1422 }

◆ Function_h()

double MonteCarloEngine::Function_h ( std::vector< double >  point)

A method to calculate the LogLikelihood + LogAprioriProbability.

Parameters
[in]pointthe set of points in the parameter space
Returns
LogLikelihood + LogAprioriProbability

Definition at line 1424 of file MonteCarloEngine.cpp.

1424  {
1425  if (point.size() != GetNParameters()) {
1426  throw std::runtime_error("MCMCENgine::Function_h : Invalid number of entries in the vector.");
1427  }
1428  return LogLikelihood(point) + LogAPrioriProbability(point);
1429 }

◆ getchainedObsSize()

int MonteCarloEngine::getchainedObsSize ( ) const
inline

A get method to access the number of Observable chains requested.

Returns
an integer for the number of Observable chains

Definition at line 195 of file MonteCarloEngine.h.

196  {
197  return kchainedObs;
198  }

◆ getHistograms1D()

std::map<std::string, BCH1D> MonteCarloEngine::getHistograms1D ( ) const
inline

A get method to access the map of 1D histograms.

Returns
the map of 1D histograms

Definition at line 242 of file MonteCarloEngine.h.

243  {
244  return Histo1D;
245  }

◆ getHistograms2D()

std::map<std::string, BCH2D> MonteCarloEngine::getHistograms2D ( ) const
inline

A get method to access the map of 2D histograms.

Returns
the map of 2D histograms

Definition at line 251 of file MonteCarloEngine.h.

252  {
253  return Histo2D;
254  }

◆ getHistoLog()

std::string MonteCarloEngine::getHistoLog ( ) const
inline

A get method to access the stream that stores the log messages coming from histogram printing and checking.

Returns
a string containing the log messages

Definition at line 186 of file MonteCarloEngine.h.

187  {
188  return HistoLog.str().c_str();
189  }

◆ getMPIWorldSize()

int MonteCarloEngine::getMPIWorldSize ( )
inline

A get method to get the number of MPI world size.

Returns
the size of the MPI world

Definition at line 422 of file MonteCarloEngine.h.

423  {
424 #ifdef _MPI
425  return procnum;
426 #else
427  return 1;
428 #endif
429  }

◆ getNumOfDiscardedEvents()

int MonteCarloEngine::getNumOfDiscardedEvents ( ) const
inline

A get method to access the number of events discarded due to failure to update model. These events are not used for the MCMC run.

Returns
the number of discarded events

Definition at line 224 of file MonteCarloEngine.h.

225  {
226  return NumOfDiscardedEvents;
227  }

◆ getNumOfUsedEvents()

int MonteCarloEngine::getNumOfUsedEvents ( ) const
inline

A get method to access the number of events used in the MCMC run.

Returns
the number of events used in the MCMC run

Definition at line 233 of file MonteCarloEngine.h.

234  {
235  return NumOfUsedEvents;
236  }

◆ getPrecision()

int MonteCarloEngine::getPrecision ( double  value,
double  rms 
)
private

A get method to calculate the precision of the numbers in the Statistics file based on the rms.

Parameters
[in]valuethe mean value
[in]rmsthe rms value
Returns
the number of digits of precision

Definition at line 1004 of file MonteCarloEngine.cpp.

1004  {
1005  if (value == 0.0) // otherwise it will return 'nan' due to the log10() of zero
1006  return 0.0;
1007 
1008  if (significants == 0) return 2 + ceil(log10(fabs(value)))-ceil(log10(rms));
1009  else return significants;
1010 }

◆ getWriteLogLikelihoodChain()

bool MonteCarloEngine::getWriteLogLikelihoodChain ( )
inline

A get method to get the value of the bool the writing of Loglikelihood in the ROOT tree.

Definition at line 333 of file MonteCarloEngine.h.

334  {
335  return WriteLogLikelihoodChain;
336  };

◆ getWriteParametersChain()

bool MonteCarloEngine::getWriteParametersChain ( )
inline

A get method to get the value of the bool the writing of parameters in the ROOT tree.

Definition at line 350 of file MonteCarloEngine.h.

351  {
352  return WriteParametersChain;
353  };

◆ InChainFillObservablesTree()

void MonteCarloEngine::InChainFillObservablesTree ( )

Definition at line 949 of file MonteCarloEngine.cpp.

950 {
951  if (!hMCMCObservableTree) return;
952  for (fMCMCTree_Chain = 0; fMCMCTree_Chain < fMCMCNChains; ++fMCMCTree_Chain) {
953  if (getchainedObsSize() > 0) hMCMCTree_Observables = hMCMCObservables[fMCMCTree_Chain];
954  if (WriteLogLikelihoodChain) {
955  hMCMCLogLikelihood = fMCMCStates.at(fMCMCTree_Chain).log_likelihood;
956  hMCMCLogProbability = fMCMCStates.at(fMCMCTree_Chain).log_probability;
957  hMCMCLogPriorProbability = fMCMCStates.at(fMCMCTree_Chain).log_prior;
958  }
959  if (kwmax > 0 && fMCMCFlagWriteChainToFile) hMCMCTree_Observables_weight = hMCMCObservables_weight[fMCMCTree_Chain];
960  hMCMCObservableTree->Fill();
961  }
962 }

◆ InChainFillParametersTree()

void MonteCarloEngine::InChainFillParametersTree ( )

Definition at line 964 of file MonteCarloEngine.cpp.

965 {
966  if (!hMCMCParameterTree) return;
967  for (fMCMCTree_Chain = 0; fMCMCTree_Chain < fMCMCNChains; ++fMCMCTree_Chain) {
968  hMCMCTree_Parameters = hMCMCParameters[fMCMCTree_Chain];
969  hMCMCParameterTree->Fill();
970  }
971 }

◆ Initialize()

void MonteCarloEngine::Initialize ( StandardModel *  Mod_i)

Initialization of the Monte Carlo Engine.

The initialization of the Monte Carlo Engine performs the following tasks

  • The observables and observables2D lists are checked and sorted into the ones which are included in the MCMC and the ones which are not.
  • The experimental likelihood for observables and observalbes2D are read if provided through a .root file.
  • The correlated Gaussian observables are sorted for inclusion in the MCMC.
    Parameters
    [in]Mod_ithe pointer to the model defined in SomeModel.conf

Definition at line 76 of file MonteCarloEngine.cpp.

77 {
78  Mod = Mod_i;
79  int k = 0, kweight = 0;
80 
81  for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin(); it < Obs_ALL.end(); it++) {
82  if (!it->isTMCMC()) {
83  k++;
84  if (it->getDistr().compare("noweight") != 0) kweight++;
85  if (it->isWriteChain()) kchainedObs++;
86  }
87  thMin[it->getName()] = std::numeric_limits<double>::max();
88  thMax[it->getName()] = -std::numeric_limits<double>::max();
89  }
90  for (std::vector<Observable2D>::iterator it = Obs2D_ALL.begin(); it < Obs2D_ALL.end(); it++) {
91  if ((it->getDistr()).compare("file") == 0) {
92  if (!it->isTMCMC())
93  throw std::runtime_error("ERROR: cannot handle noMCMC for Observable2D file yet!");
94  } else if (it->getDistr().compare("weight") == 0)
95  throw std::runtime_error("ERROR: do not use Observable2D for analytic 2D weights!");
96  }
97  for (std::vector<CorrelatedGaussianObservables>::iterator it1 = CGO.begin(); it1 != CGO.end(); ++it1) {
98  std::vector<Observable> ObsV(it1->getObs());
99  for (std::vector<Observable>::iterator it = ObsV.begin(); it != ObsV.end(); ++it) {
100  if ((it->getDistr()).compare("file") == 0)
101  throw std::runtime_error("Cannot use file in CorrelatedGaussianObservables!");
102  if (!(it->isTMCMC())) {
103  k++;
104  if (it->getDistr().compare("noweight") != 0)
105  throw std::runtime_error("Cannot use weight in CorrelatedGaussianObservables!");
106  }
107  thMin[it->getName()] = std::numeric_limits<double>::max();
108  thMax[it->getName()] = -std::numeric_limits<double>::max();
109  }
110  if (it1->isPrediction()) {
111  CorrelationMap[it1->getName()] = new TPrincipal(it1->getObs().size(), "N");
112  }
113  }
114  kmax = k;
115  kwmax = kweight;
116 
117  unknownParameters = Mod->getUnknownParameters();
118  DefineParameters();
119 };

◆ LogLikelihood()

double MonteCarloEngine::LogLikelihood ( const std::vector< double > &  parameters)

This member calculates the loglikelihood for the observables included in the MCMC run.

The model is updated with the new set of parameters through the Model::Update() method. If this is successful this even is counted as a used event. Otherwise it is counted as a discarded event. The log probability is calculated using Weight(const Observable&, const double& ) function for the observables, Weight(const Observable2D&, const double&, const double& ) for the observable2D and Weight(const CorrelatedGaussianObservables& ) for the correlated Gaussian observables. Overloaded from BCEngineMCMC in BAT.

Parameters
[in]parametersthe vector of the parameters that are being varied in the Monte Carlo
Returns
the loglikelihood

Definition at line 312 of file MonteCarloEngine.cpp.

312  {
313  // This methods returns the logarithm of the conditional probability
314  // p(data|parameters). This is where you have to define your model.
315 
316  double logprob = 0.;
317 
318  setDParsFromParameters(parameters, DPars);
319 
320  // if update false set probability equal zero
321  if (!Mod->Update(DPars)) {
322 #ifdef _MCDEBUG
323  std::cout << "event discarded" << std::endl;
324 
325  /* Debug */
326  //for (int k = 0; k < parameters.size(); k++)
327  // std::cout << " " << GetParameter(k)->GetName() << " = "
328  // << DPars[GetParameter(k)->GetName()] << std::endl;
329 #endif
330  NumOfDiscardedEvents++;
331  return (log(0.));
332  }
333 #ifdef _MCDEBUG
334  //std::cout << "event used in MC" << std::endl;
335 #endif
336 
337  for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin(); it != Obs_ALL.end(); it++) {
338  if (it->isTMCMC()) logprob += it->computeWeight();
339  }
340 
341  for (std::vector<Observable2D>::iterator it = Obs2D_ALL.begin(); it != Obs2D_ALL.end(); it++) {
342  if (it->isTMCMC()) logprob += it->computeWeight();
343  }
344 
345  for (std::vector<CorrelatedGaussianObservables>::iterator it = CGO.begin(); it < CGO.end(); it++) {
346  if(!(it->isPrediction())) logprob += it->computeWeight();
347  }
348  if (!std::isfinite(logprob)) {
349  NumOfDiscardedEvents++;
350 #ifdef _MCDEBUG
351 // std::cout << "Event discarded since logprob evaluated to: " << logprob << std::endl ;
352 #endif
353  return (log(0.));
354  }
355  NumOfUsedEvents++;
356  return logprob;
357 }

◆ MCMCUserIterationInterface()

void MonteCarloEngine::MCMCUserIterationInterface ( )

Overloaded from BCEngineMCMC in BAT

The interface is used to update the model parameters using the Model::Update() method. Then the values of the observables are computed and the respective histograms are filled. If the checkTheoryRange flag is true then the minimum and maximum of the theory value is checked and reset to include the current theory values.

Definition at line 359 of file MonteCarloEngine.cpp.

359  {
360 #ifdef _MPI
361  unsigned mychain = 0;
362  int iproc = 0;
363  unsigned npars = GetNParameters();
364  int buffsize = npars + 1;
365  int index_chain[procnum];
366  double *recvbuff = new double[buffsize];
367  std::vector<double> pars;
368  double **buff;
369 
370  buff = new double*[procnum];
371  int obsbuffsize = Obs_ALL.size() + 2 * Obs2D_ALL.size();
372  for (std::vector<CorrelatedGaussianObservables>::iterator it1 = CGO.begin(); it1 < CGO.end(); it1++)
373  obsbuffsize += it1->getObs().size();
374  buff[0] = new double[procnum * obsbuffsize];
375  for (int i = 1; i < procnum; i++) {
376  buff[i] = buff[i - 1] + obsbuffsize;
377  index_chain[i] = -1;
378  }
379 
380  double ** sendbuff = new double *[procnum];
381  sendbuff[0] = new double[procnum * buffsize];
382  for (int il = 1; il < procnum; il++)
383  sendbuff[il] = sendbuff[il - 1] + buffsize;
384 
385  while (mychain < fMCMCNChains) {
386  pars.clear();
387  pars = fMCMCStates.at(mychain).parameters;
388 
389  if (PrintLoglikelihoodPlots || WriteParametersChain) {
390  std::map<std::string, double> tmpDPars;
391  setDParsFromParameters(pars, tmpDPars);
392  if (PrintLoglikelihoodPlots) DPars_allChains.push_back(tmpDPars);
393  if (WriteParametersChain) {
394  int k = 0;
395  for (std::map<std::string, double>::iterator it = tmpDPars.begin(); it != tmpDPars.end(); it++) hMCMCParameters[mychain][k++] = it->second;
396  }
397  }
398 
399  index_chain[iproc] = mychain;
400  iproc++;
401  mychain++;
402  if (iproc < procnum && mychain < fMCMCNChains)
403  continue;
404 
405  for (int il = 0; il < iproc; il++) {
406  //The first entry of the array specifies the task to be executed.
407 
408  sendbuff[il][0] = 2.; // 2 = observables calculation
409  for (int im = 1; im < buffsize; im++) sendbuff[il][im] = fMCMCStates.at(index_chain[il]).parameters.at(im - 1);
410  }
411  for (int il = iproc; il < procnum; il++) {
412  sendbuff[il][0] = 3.; // 3 = nothing to execute, but return a buffer of observables
413  index_chain[il] = -1;
414  }
415  // double inittime = MPI::Wtime();
416  MPI_Scatter(sendbuff[0], buffsize, MPI_DOUBLE,
417  recvbuff, buffsize, MPI_DOUBLE,
418  0, MPI_COMM_WORLD);
419 
420  if (recvbuff[0] == 2.) { // compute observables
421  double sbuff[obsbuffsize];
422  pars.assign(recvbuff + 1, recvbuff + buffsize);
423  setDParsFromParameters(pars,DPars);
424  Mod->Update(DPars);
425 
426  int k = 0;
427  for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin(); it < Obs_ALL.end(); it++) {
428  sbuff[k++] = it->computeTheoryValue();
429  }
430  for (std::vector<Observable2D>::iterator it = Obs2D_ALL.begin(); it < Obs2D_ALL.end(); it++) {
431  sbuff[k++] = it->computeTheoryValue();
432  sbuff[k++] = it->computeTheoryValue2();
433  }
434 
435  for (std::vector<CorrelatedGaussianObservables>::iterator it1 = CGO.begin(); it1 < CGO.end(); it1++) {
436  std::vector<Observable> ObsV(it1->getObs());
437  for (std::vector<Observable>::iterator it = ObsV.begin(); it != ObsV.end(); ++it)
438  sbuff[k++] = it->computeTheoryValue();
439  }
440  MPI_Gather(sbuff, obsbuffsize, MPI_DOUBLE, buff[0], obsbuffsize, MPI_DOUBLE, 0, MPI_COMM_WORLD);
441  } else if (recvbuff[0] == 3.) { // do not compute observables, but gather the buffer
442  double sbuff[obsbuffsize];
443  MPI_Gather(sbuff, obsbuffsize, MPI_DOUBLE, buff[0], obsbuffsize, MPI_DOUBLE, 0, MPI_COMM_WORLD);
444  }
445 
446  for (int il = 0; il < procnum; il++) {
447  if (index_chain[il] >= 0) {
448  int k = 0;
449  // fill the histograms for observables
450  int ko = 0, kweight = 0, k_all = 0, k_cObs = 0;
451  for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin();
452  it < Obs_ALL.end(); it++) {
453  double th = buff[il][k++];
454  /* set the min and max of theory values */
455  if (th < thMin[it->getName()]) thMin[it->getName()] = th;
456  if (th > thMax[it->getName()]) thMax[it->getName()] = th;
457  Histo1D[it->getName()].GetHistogram()->Fill(th);
458  if (fMCMCFlagWriteChainToFile) hMCMCObservables[index_chain[il]][k_all++] = th;
459  else if (getchainedObsSize() > 0 && it->isWriteChain()) hMCMCObservables[index_chain[il]][k_cObs++] = th;
460  if (!it->isTMCMC()) {
461  obval[index_chain[il] * kmax + ko] = th;
462  ko++;
463  if (it->getDistr().compare("noweight") != 0 && it->getDistr().compare("writeChain") != 0) {
464  double weight = it->computeWeight(th);
465  obweight[index_chain[il] * kwmax + kweight] = weight;
466  if (fMCMCFlagWriteChainToFile) hMCMCObservables_weight[index_chain[il]][kweight] = weight;
467  kweight++;
468  }
469  }
470  }
471 
472  // fill the 2D histograms for observables
473  for (std::vector<Observable2D>::iterator it = Obs2D_ALL.begin();
474  it < Obs2D_ALL.end(); it++) {
475  double th1 = buff[il][k++];
476  double th2 = buff[il][k++];
477  Histo2D[it->getName()].GetHistogram()->Fill(th1, th2);
478  }
479 
480  // fill the histograms for correlated observables
481  for (std::vector<CorrelatedGaussianObservables>::iterator it1 = CGO.begin(); it1 < CGO.end(); it1++) {
482  std::vector<Observable> ObsV(it1->getObs());
483  Double_t * COdata = new Double_t[ObsV.size()];
484  int nObs = 0;
485  for (std::vector<Observable>::iterator it = ObsV.begin(); it != ObsV.end(); ++it) {
486  double th = buff[il][k++];
487  /* set the min and max of theory values */
488  if (th < thMin[it->getName()]) thMin[it->getName()] = th;
489  if (th > thMax[it->getName()]) thMax[it->getName()] = th;
490  Histo1D[it->getName()].GetHistogram()->Fill(th);
491  if (it1->isPrediction()) COdata[nObs++] = th;
492  }
493  if (it1->isPrediction()) CorrelationMap[it1->getName()]->AddRow(COdata);
494  delete [] COdata;
495  }
496  }
497  }
498  iproc = 0;
499  }
500  if (fMCMCFlagWriteChainToFile || getchainedObsSize() > 0 || WriteLogLikelihoodChain) InChainFillObservablesTree();
501  if (WriteParametersChain) InChainFillParametersTree();
502  delete sendbuff[0];
503  delete [] sendbuff;
504  delete [] recvbuff;
505  delete buff[0];
506  delete [] buff;
507 #else
508  for (unsigned int i = 0; i < fMCMCNChains; ++i) {
509  // NOTE: BAT syncs fMCMCThreadLocalStorage with fMCMCStates before calling MCMCUserIterationInterface.
510  std::vector<double>::const_iterator first = fMCMCStates.at(i).parameters.begin();
511  std::vector<double>::const_iterator last = first + GetNParameters();
512  std::vector<double> currvec(first, last);
513  setDParsFromParameters(currvec,DPars);
514  if (PrintLoglikelihoodPlots) DPars_allChains.push_back(DPars);
515  if (WriteParametersChain) {
516  int k = 0;
517  for (std::map<std::string, double>::iterator it = DPars.begin(); it != DPars.end(); it++) hMCMCParameters[i][k++] = it->second;
518  }
519 
520  Mod->Update(DPars);
521  // fill the histograms for observables
522  int k = 0, kweight = 0, k_all = 0, k_cObs = 0;
523  for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin();
524  it < Obs_ALL.end(); it++) {
525  double th = it->computeTheoryValue();
526  /* set the min and max of theory values */
527  if (th < thMin[it->getName()]) thMin[it->getName()] = th;
528  if (th > thMax[it->getName()]) thMax[it->getName()] = th;
529  Histo1D[it->getName()].GetHistogram()->Fill(th);
530  if (fMCMCFlagWriteChainToFile) hMCMCObservables[i][k_all++] = th;
531  else if (getchainedObsSize() > 0 && it->isWriteChain()) hMCMCObservables[i][k_cObs++] = th;
532  if (!it->isTMCMC()) {
533  obval[i * kmax + k] = th;
534  k++;
535  if (it->getDistr().compare("noweight") != 0 && it->getDistr().compare("writeChain") != 0) {
536  double weight = it->computeWeight(th);
537  obweight[i * kwmax + kweight] = weight;
538  if (fMCMCFlagWriteChainToFile) hMCMCObservables_weight[i][kweight] = weight;
539  kweight++;
540  }
541  }
542  }
543 
544  // fill the 2D histograms for observables
545  for (std::vector<Observable2D>::iterator it = Obs2D_ALL.begin();
546  it < Obs2D_ALL.end(); it++) {
547  double th1 = it->computeTheoryValue();
548  double th2 = it->computeTheoryValue2();
549  Histo2D[it->getName()].GetHistogram()->Fill(th1, th2);
550  }
551 
552  // fill the histograms for correlated observables
553  for (std::vector<CorrelatedGaussianObservables>::iterator it1 = CGO.begin();
554  it1 < CGO.end(); it1++) {
555  std::vector<Observable> ObsV(it1->getObs());
556  Double_t * COdata = new Double_t[ObsV.size()];
557  int nObs = 0;
558  for (std::vector<Observable>::iterator it = ObsV.begin();
559  it != ObsV.end(); ++it) {
560  double th = it->computeTheoryValue();
561  /* set the min and max of theory values */
562  if (th < thMin[it->getName()]) thMin[it->getName()] = th;
563  if (th > thMax[it->getName()]) thMax[it->getName()] = th;
564  Histo1D[it->getName()].GetHistogram()->Fill(th);
565  if (it1->isPrediction()) COdata[nObs++] = th;
566  }
567  if (it1->isPrediction()) CorrelationMap[it1->getName()]->AddRow(COdata);
568  delete [] COdata;
569  }
570  }
571 
572  if (fMCMCFlagWriteChainToFile || getchainedObsSize() > 0 || WriteLogLikelihoodChain) InChainFillObservablesTree();
573  if (WriteParametersChain) InChainFillParametersTree();
574 #endif
575  for (unsigned int i = 0; i < fMCMCNChains; i++) {
576  double LogLikelihood = fMCMCStates.at(i).log_likelihood;
577  if (LogLikelihood > LogLikelihood_max) {
578  LogLikelihood_max = std::max(LogLikelihood, LogLikelihood_max);
579  par_at_LL_max = Getx(i); // NOTE: Unrotated, to be rotated later.
580  }
581  Histo1D["LogLikelihood"].GetHistogram()->Fill(LogLikelihood);
582  if (PrintLoglikelihoodPlots) {
583  for (std::vector<ModelParameter>::const_iterator it = ModPars.begin(); it != ModPars.end(); it++) {
584  if (it->IsFixed()) continue;
585  if (std::find(unknownParameters.begin(), unknownParameters.end(), it->getname()) != unknownParameters.end()) continue;
586  std::string HistName = it->getname() + "_vs_LogLikelihood";
587  Histo2D[HistName].GetHistogram()->Fill(DPars_allChains.at(i).at(it->getname()), LogLikelihood);
588  }
589  }
590  }
591  if (PrintLoglikelihoodPlots) DPars_allChains.clear();
592 }

◆ Print1D()

void MonteCarloEngine::Print1D ( BCH1D  hist1D,
const char *  filename,
int  ww = 0,
int  wh = 0 
)

Definition at line 619 of file MonteCarloEngine.cpp.

619  {
620  TCanvas * c;
621  cindex++;
622  if(ww > 0 && wh > 0)
623  c = new TCanvas(TString::Format("c_bch1d_%d",cindex), TString::Format("c_bch1d_%d",cindex), ww, wh);
624  else
625  c = new TCanvas(TString::Format("c_bch1d_%d",cindex));
626 
627  bch1d.GetHistogram()->Scale(1./bch1d.GetHistogram()->Integral("width"));
628 
629  bch1d.SetBandType(BCH1D::kSmallestInterval);
630  bch1d.SetBandColor(0, gIdx);
631  bch1d.SetBandColor(1, rIdx);
632  bch1d.SetBandColor(2, kOrange - 3);
633  bch1d.SetNBands(3);
634  bch1d.SetNSmooth(nSmooth);
635  bch1d.SetDrawGlobalMode(true);
636  bch1d.SetDrawMean(true, true);
637  bch1d.SetDrawLegend(!noLegend);
638  if (noLegend) gStyle->SetOptStat("emr");
639  bch1d.SetNLegendColumns(1);
640  bch1d.SetStats(true);
641 
642  bch1d.Draw();
643 
644  if (printLogo) {
645  double xRange = (bch1d.GetHistogram()->GetXaxis()->GetXmax() - bch1d.GetHistogram()->GetXaxis()->GetXmin())*3./4.;
646  double yRange = (bch1d.GetHistogram()->GetMaximum() - bch1d.GetHistogram()->GetMinimum());
647 
648  double xL;
649  if (noLegend) xL = bch1d.GetHistogram()->GetXaxis()->GetXmin()+0.0475*xRange;
650  else xL = bch1d.GetHistogram()->GetXaxis()->GetXmin() + 0.0375 * xRange;
651  double yL = bch1d.GetHistogram()->GetYaxis()->GetXmin() + 0.89 * yRange;
652 
653  double xR;
654  if (noLegend) xR = xL + 0.21*xRange;
655  else xR = xL + 0.18 * xRange;
656  double yR = yL + 0.09 * yRange;
657 
658  TBox b1 = TBox(xL, yL, xR, yR);
659  b1.SetFillColor(gIdx);
660 
661  TBox b2;
662  b2 = TBox(xL+0.008*xRange, yL+0.008*yRange, xR-0.008*xRange, yR-0.008*yRange);
663  b2.SetFillColor(kWhite);
664 
665  TPaveText b3 = TPaveText(xL+0.014*xRange, yL+0.013*yRange, xL+0.70*(xR-xL), yR-0.013*yRange);
666  if (noLegend) b3.SetTextSize(0.056);
667  else b3.SetTextSize(0.051);
668  b3.SetTextAlign(22);
669  b3.SetTextColor(kWhite);
670  b3.AddText("HEP");
671  b3.SetFillColor(rIdx);
672 
673  TPaveText * b4;
674  if (noLegend) {
675  b4 = new TPaveText(xL+0.72*(xR-xL), yL+0.030*yRange, xR-0.008*xRange, yR-0.013*yRange);
676  b4->SetTextSize(0.048);
677  } else {
678  b4 = new TPaveText(xL + 0.75 * (xR - xL), yL + 0.024 * yRange, xR - 0.008 * xRange, yR - 0.013 * yRange);
679  b4->SetTextSize(0.039);
680  }
681  b4->SetTextAlign(33);
682  b4->SetTextColor(rIdx);
683  b4->AddText("fit");
684  b4->SetFillColor(kWhite);
685 
686  b1.Draw("SAME");
687  b2.Draw("SAME");
688  b3.Draw("SAME");
689  b4->Draw("SAME");
690 
691  c->Print(filename);
692  delete b4;
693  b4 = NULL;
694  } else c->Print(filename);
695 
696  delete c;
697  c = NULL;
698 }

◆ Print2D()

void MonteCarloEngine::Print2D ( BCH2D  hist2D,
const char *  filename,
int  ww = 0,
int  wh = 0 
)

Definition at line 700 of file MonteCarloEngine.cpp.

701 {
702  TCanvas * c;
703  cindex++;
704  if(ww > 0 && wh > 0)
705  c = new TCanvas(TString::Format("c_bch2d_%d",cindex), TString::Format("c_bch2d_%d",cindex), ww, wh);
706  else
707  c = new TCanvas(TString::Format("c_bch2d_%d",cindex));
708 
709  bch2d.GetHistogram()->Scale(1./bch2d.GetHistogram()->Integral("width"));
710  bch2d.GetHistogram()->GetYaxis()->SetTitleOffset(1.45);
711 
712  bch2d.SetBandType(BCH2D::kSmallestInterval);
713  bch2d.SetBandColor(0, TColor::GetColorTransparent(kOrange - 3, alpha2D));
714  bch2d.SetBandColor(1, TColor::GetColorTransparent(rIdx, alpha2D));
715  bch2d.SetBandColor(2, TColor::GetColorTransparent(gIdx, alpha2D));
716  bch2d.SetNBands(3);
717  bch2d.SetBandFillStyle(histogram2Dtype);// Type of 2D Histogram 1001 -> box pixel, 101 -> filled, 1 -> contour.
718  if (histogram2Dtype == 1 || histogram2Dtype == 101) bch2d.SetNSmooth(1);
719  else bch2d.SetNSmooth(0);
720  if (histogram2Dtype == 1) bch2d.GetHistogram()->SetLineWidth(3);
721  bch2d.SetDrawLocalMode(false);
722  bch2d.SetDrawGlobalMode(true);
723  bch2d.SetDrawMean(true, true);
724  bch2d.SetDrawLegend(!noLegend);
725  if (noLegend) gStyle->SetOptStat("emr");
726  bch2d.SetNLegendColumns(1);
727  bch2d.SetStats(true);
728 
729  bch2d.Draw();
730 
731  if (printLogo) {
732  double xRange = (bch2d.GetHistogram()->GetXaxis()->GetXmax() - bch2d.GetHistogram()->GetXaxis()->GetXmin())*3./4.;
733  double yRange = bch2d.GetHistogram()->GetYaxis()->GetXmax() - bch2d.GetHistogram()->GetYaxis()->GetXmin();
734 
735  double xL;
736  if (noLegend) xL = bch2d.GetHistogram()->GetXaxis()->GetXmin()+0.0475*xRange;
737  else xL = bch2d.GetHistogram()->GetXaxis()->GetXmin() + 0.0375 * xRange;
738  double yL = bch2d.GetHistogram()->GetYaxis()->GetXmin() + 0.89 * yRange;
739 
740  double xR;
741  if (noLegend) xR = xL + 0.21*xRange;
742  else xR = xL + 0.18 * xRange;
743  double yR = yL + 0.09 * yRange;
744 
745  TBox b1 = TBox(xL, yL, xR, yR);
746  b1.SetFillColor(gIdx);
747 
748  TBox b2 = TBox(xL + 0.008 * xRange, yL + 0.008 * yRange, xR - 0.008 * xRange, yR - 0.008 * yRange);
749  b2.SetFillColor(kWhite);
750 
751  TPaveText b3 = TPaveText(xL + 0.014 * xRange, yL + 0.013 * yRange, xL + 0.70 * (xR - xL), yR - 0.013 * yRange);
752  if (noLegend) b3.SetTextSize(0.056);
753  else b3.SetTextSize(0.051);
754  b3.SetTextAlign(22);
755  b3.SetTextColor(kWhite);
756  b3.AddText("HEP");
757  b3.SetFillColor(rIdx);
758 
759  TPaveText * b4;
760  if (noLegend) {
761  b4 = new TPaveText(xL+0.72*(xR-xL), yL+0.030*yRange, xR-0.008*xRange, yR-0.013*yRange);
762  b4->SetTextSize(0.048);
763  } else {
764  b4 = new TPaveText(xL + 0.75 * (xR - xL), yL + 0.024 * yRange, xR - 0.008 * xRange, yR - 0.013 * yRange);
765  b4->SetTextSize(0.039);
766  }
767  b4->SetTextAlign(33);
768  b4->SetTextColor(rIdx);
769  b4->AddText("fit");
770  b4->SetFillColor(kWhite);
771 
772  b1.Draw("SAME");
773  b2.Draw("SAME");
774  b3.Draw("SAME");
775  b4->Draw("SAME");
776 
777  c->Print(filename);
778  delete b4;
779  b4 = NULL;
780  } else c->Print(filename);
781 
782  delete c;
783  c = NULL;
784 }

◆ PrintCorrelationMatrixToLaTeX()

void MonteCarloEngine::PrintCorrelationMatrixToLaTeX ( const std::string  filename)

This member generates the correlation matrix using BCH2D from the BAT libraries.

Parameters
[in]filenamethe name of the file where the correlation matrix is printed

Definition at line 973 of file MonteCarloEngine.cpp.

973  {
974  std::ofstream out;
975  out.open(filename.c_str(), std::ios::out);
976 
977  int npar = GetNParameters();
978 
979  for (int i = 0; i < npar; ++i)
980  out << " & " << GetParameter(i).GetName();
981  out << " \\\\" << std::endl;
982 
983  for (int i = 0; i < npar; ++i) {
984  out << GetParameter(i).GetName() << " & $";
985  for (int j = 0; j < npar; ++j) {
986  if (i != j) {
987  BCH2D* bch2d_temp = new BCH2D(GetMarginalized(GetParameter(i).GetName(), GetParameter(j).GetName()));
988  if (bch2d_temp != NULL)
989  out << bch2d_temp->GetHistogram()->GetCorrelationFactor();
990  else
991  out << 0.;
992  delete bch2d_temp;
993  bch2d_temp = NULL;
994  } else
995  out << 1.;
996  if (j == npar - 1) out << "$ \\\\" << std::endl;
997  else out << "$ & $";
998  }
999  }
1000 
1001  out.close();
1002 }

◆ PrintHistogram() [1/2]

void MonteCarloEngine::PrintHistogram ( std::string &  OutFile,
const std::string  OutputDir 
)

Member used for printing histograms for observables, observable2D, correlated Gaussian observables and model parameters vs. observables.

Parameters
[in]outa reference to an object of type BCModelOutput as defined in the BAT libraries
[in]OutputDirthe name of the output directory

Definition at line 814 of file MonteCarloEngine.cpp.

815 {
816  std::vector<double> mode(GetBestFitParameters());
817  if (mode.size() == 0) throw std::runtime_error("\n ERROR: Global Mode could not be determined possibly because of infinite loglikelihood. Observables histogram cannot be generated.\n");
818  setDParsFromParameters(mode,DPars);
819 
820  Mod->Update(DPars);
821 
822  if (Obs_ALL.size() != 0 || CGO.size() != 0) std::cout << "\nPrinting 1D histograms in the Observables directory: " << std::endl;
823  for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin(); it < Obs_ALL.end(); it++) PrintHistogram(OutFile, *it, OutputDir);
824 
825  for (std::vector<CorrelatedGaussianObservables>::iterator it1 = CGO.begin(); it1 < CGO.end(); it1++) {
826  std::vector<Observable> ObsV(it1->getObs());
827  for (std::vector<Observable>::iterator it = ObsV.begin(); it != ObsV.end(); ++it) PrintHistogram(OutFile, *it, OutputDir);
828  }
829 
830  if (Obs2D_ALL.size() != 0) std::cout << "\nPrinting 2D histograms in the Observables directory: " << std::endl;
831  for (std::vector<Observable2D>::iterator it = Obs2D_ALL.begin(); it < Obs2D_ALL.end(); it++) {
832  std::string HistName = it->getName();
833  if (Histo2D[HistName].GetHistogram()->Integral() > 0.0) {
834  std::string fname = OutputDir + "/" + HistName + ".pdf";
835  std::vector<double> th;
836  th.push_back(it->computeTheoryValue());
837  th.push_back(it->computeTheoryValue2());
838  Histo2D[HistName].SetGlobalMode(th);
839  Print2D(Histo2D[HistName], fname.c_str());
840  std::cout << " " + HistName + ".pdf" << std::endl;
841  if(OutFile.compare("") == 0) throw std::runtime_error("\nMonteCarloEngine::PrintHistogram ERROR: No root file specified for writing histograms.");
842  TDirectory * dir = gDirectory;
843  GetOutputFile()->cd();
844  Histo2D[HistName].GetHistogram()->Write();
845  gDirectory = dir;
846  CheckHistogram(*Histo2D[HistName].GetHistogram(), HistName);
847  } else HistoLog << "WARNING: The histogram of " << HistName << " is empty!" << std::endl;
848  }
849 
850  std::cout << "\nPrinting LogLikelihood histogram in the Observables directory: " << std::endl;
851  if (Histo1D["LogLikelihood"].GetHistogram()->Integral() > 0.0) {
852  std::string fname = OutputDir + "/LogLikelihood.pdf";
853  Print1D(Histo1D["LogLikelihood"], fname.c_str());
854  std::cout << " LogLikelihood.pdf" << std::endl;
855  TDirectory * dir = gDirectory;
856  GetOutputFile()->cd();
857  Histo1D["LogLikelihood"].GetHistogram()->Write();
858  gDirectory = dir;
859  CheckHistogram(*Histo1D["LogLikelihood"].GetHistogram(), "LogLikelihood");
860  }
861 
862  if (PrintLoglikelihoodPlots) {
863  std::cout << "\nPrinting LogLikelihood vs. parameter 2D histograms in the Observables directory: " << std::endl;
864  for (std::vector<ModelParameter>::const_iterator it = ModPars.begin(); it != ModPars.end(); it++) {
865  if (it->IsFixed()) continue;
866  if (std::find(unknownParameters.begin(), unknownParameters.end(), it->getname()) != unknownParameters.end()) continue;
867  std::string HistName = it->getname() + "_vs_LogLikelihood";
868  if (Histo2D[HistName].GetHistogram()->Integral() > 0.0) {
869  std::string fname = OutputDir + "/LogLikelihoodPlots/" + HistName + ".pdf";
870  Print2D(Histo2D[HistName], fname.c_str());
871  std::cout << " " + HistName + ".pdf" << std::endl;
872  if (OutFile.compare("") == 0) throw std::runtime_error("\nMonteCarloEngine::PrintHistogram ERROR: No root file specified for writing histograms.");
873  TDirectory * dir = gDirectory;
874  GetOutputFile()->cd();
875  Histo2D[HistName].GetHistogram()->Write();
876  gDirectory = dir;
877  CheckHistogram(*Histo2D[HistName].GetHistogram(), HistName);
878  } else HistoLog << "WARNING: The histogram of " << HistName << " is empty!" << std::endl;
879  }
880  }
881  if (noLegend) gStyle->SetOptStat("emrn");
882 }

◆ PrintHistogram() [2/2]

void MonteCarloEngine::PrintHistogram ( std::string &  OutFile,
Observable &  it,
const std::string  OutputDir 
)

Overloaded from PrintHistogram(BCModelOutput&, const std::string) to print histogram for observables.

Parameters
[in]outa reference to an object of type BCModelOutput as defined in the BAT libraries
[in]ita reference to an object of type Observable
[in]OutputDirthe name of the output directory

Definition at line 786 of file MonteCarloEngine.cpp.

787 {
788 
789  std::string HistName = it.getName();
790  double min = thMin[it.getName()];
791  double max = thMax[it.getName()];
792  if (Histo1D[HistName].GetHistogram()->Integral() > 0.0) {
793  std::string fname = OutputDir + "/" + HistName + ".pdf";
794  Histo1D[HistName].SetGlobalMode(it.computeTheoryValue());
795  Print1D(Histo1D[HistName], fname.c_str());
796  std::cout << " " + HistName + ".pdf" << std::endl;
797  if(OutFile.compare("") == 0) {
798  throw std::runtime_error("\nMonteCarloEngine::PrintHistogram ERROR: No root file specified for writing histograms.");
799  }
800  TDirectory * dir = gDirectory;
801  GetOutputFile()->cd();
802  Histo1D[HistName].GetHistogram()->Write();
803  gDirectory = dir;
804  CheckHistogram(*Histo1D[HistName].GetHistogram(), it.getName());
805  } else
806  HistoLog << "WARNING: The histogram of "
807  << it.getName() << " is empty!" << std::endl;
808 
809  HistoLog.precision(10);
810  HistoLog << " [min, max]=[" << min << ", " << max << "]" << std::endl;
811  HistoLog.precision(6);
812 }

◆ SecondDerivative()

double MonteCarloEngine::SecondDerivative ( BCParameter  par1,
BCParameter  par2,
std::vector< double >  point 
)

A method to calculate the second derivative.

Returns
the second derivative

Definition at line 1354 of file MonteCarloEngine.cpp.

1354  {
1355 
1356  if (point.size() != GetNParameters()) {
1357  throw std::runtime_error("MCMCENgine::SecondDerivative : Invalid number of entries in the vector.");
1358  }
1359 
1360  // define steps
1361  const double dy1 = par2.GetRangeWidth() / NSTEPS;
1362  const double dy2 = dy1 * 2.;
1363  const double dy3 = dy1 * 3.;
1364 
1365  // define points at which to evaluate
1366  std::vector<double> y1p = point;
1367  std::vector<double> y1m = point;
1368  std::vector<double> y2p = point;
1369  std::vector<double> y2m = point;
1370  std::vector<double> y3p = point;
1371  std::vector<double> y3m = point;
1372 
1373  unsigned idy = GetParameters().Index(par2.GetName());
1374 
1375  y1p[idy] += dy1;
1376  y1m[idy] -= dy1;
1377  y2p[idy] += dy2;
1378  y2m[idy] -= dy2;
1379  y3p[idy] += dy3;
1380  y3m[idy] -= dy3;
1381 
1382  const double m1 = (FirstDerivative(par1, y1p) - FirstDerivative(par1, y1m)) / 2. / dy1;
1383  const double m2 = (FirstDerivative(par1, y2p) - FirstDerivative(par1, y2m)) / 4. / dy1;
1384  const double m3 = (FirstDerivative(par1, y3p) - FirstDerivative(par1, y3m)) / 6. / dy1;
1385 
1386  return 3. / 2. * m1 - 3. / 5. * m2 + 1. / 10. * m3;
1387 }

◆ setAlpha2D()

void MonteCarloEngine::setAlpha2D ( double  alpha)
inline

A set method to toggle the printing of legends in 1D and 2D histograms.

Parameters
[in]legenda boolean to toggle the the printing of legends

Definition at line 359 of file MonteCarloEngine.h.

360  {
361  alpha2D = alpha;
362  };

◆ setDParsFromParameters()

void MonteCarloEngine::setDParsFromParameters ( const std::vector< double > &  parameters,
std::map< std::string, double > &  DPars_i 
)

A method to rotate the diagonalized parameters to the original basis for correlated parameters.

Parameters
[in]parametersthe parameters in the diagonal basis
[in]DPars_ithe parameters in the in the input basis

Definition at line 261 of file MonteCarloEngine.cpp.

263 {
264  std::map<std::string, std::vector<double> > cgpmap;
265 
266  unsigned int k = 0;
267  for (std::vector<ModelParameter>::const_iterator it = ModPars.begin(); it != ModPars.end(); it++){
268  if(it->IsFixed())
269  continue;
270  if (std::find(unknownParameters.begin(), unknownParameters.end(), it->getname()) != unknownParameters.end())
271  continue;
272  if(it->getname().compare(GetParameter(k).GetName()) != 0)
273  {
274  std::stringstream out;
275  out << it->getname();
276  throw std::runtime_error("MonteCarloEngine::setDParsFromParameters(): " + out.str() + "is sitting at the wrong position in the BAT parameters vector");
277  }
278  if (it->IsCorrelated()) {
279  std::string index = it->getname().substr(it->getCgp_name().size());
280  unsigned long int lindex = strtol(index.c_str(),NULL,10);
281  if (lindex - 1 == cgpmap[it->getCgp_name()].size())
282  cgpmap[it->getCgp_name()].push_back(parameters[k]);
283  else {
284  std::stringstream out;
285  out << it->getname() << " " << lindex;
286  throw std::runtime_error("MonteCarloEngine::setDParsFromParameters(): " + out.str() + "seems to be a CorrelatedGaussianParameters object but the corresponding parameters are missing or not in the right order");
287  }
288 
289  } else
290  DPars_i[it->getname()] = parameters[k];
291  k++;
292  }
293 
294  for (unsigned int j = 0; j < CGP.size(); j++) {
295  std::vector<double> current = cgpmap.at(CGP[j].getName());
296  if (current.size() != CGP[j].getPars().size()) {
297  std::stringstream out;
298  out << CGP[j].getName();
299  throw std::runtime_error("MonteCarloEngine::setDParsFromParameters(): " + out.str() + " appears to be represented in cgpmap with a wrong size");
300  }
301 
302  std::vector<double> porig = CGP[j].getOrigParsValue(current);
303 
304  for(unsigned int l = 0; l < porig.size(); l++) {
305  DPars_i[CGP[j].getPar(l).getname()] = porig[l];
306  }
307  }
308 }

◆ setHistogram2DType()

void MonteCarloEngine::setHistogram2DType ( int  type)
inline

A set method to set the band fill type for 2D histograms.

Parameters
[in]type2D Histogram type, 1001 -> Lego, 101 -> Filled, 0 -> Contour

Definition at line 377 of file MonteCarloEngine.h.

378  {
379  histogram2Dtype = type;
380  };

◆ setHistogramBufferSize()

void MonteCarloEngine::setHistogramBufferSize ( unsigned int  histogramBufferSize_i)
inline

A set method to set the size of the buffer used by the histograms.

Parameters
[in]histogramBufferSize_ithe size of the buffer used by the histograms

Definition at line 413 of file MonteCarloEngine.h.

414  {
415  histogramBufferSize = histogramBufferSize_i;
416  };

◆ setNBins1D()

void MonteCarloEngine::setNBins1D ( unsigned int  nbins)
inline

A set method to set the number of bins for a 1D histograms.

Parameters
[in]nbinsthe number of bins in the 1D histogram

Definition at line 386 of file MonteCarloEngine.h.

387  {
388  nBins1D = nbins;
389  };

◆ setNBins2D()

void MonteCarloEngine::setNBins2D ( unsigned int  nbins)
inline

A set method to set the number of bins for a 2D histograms.

Parameters
[in]nbinsthe number of bins in the 2D histogram

Definition at line 395 of file MonteCarloEngine.h.

396  {
397  nBins2D = nbins;
398  };

◆ setNChains()

void MonteCarloEngine::setNChains ( unsigned int  i)

A set method to fix the number of chains.

The number of chains are set using the MCMCSetNChains() from BCEngineMCMC class in BAT. This also creates pointers to the vector of observable values and weights for all the chains.

Parameters
[in]ithe number of chains

Definition at line 178 of file MonteCarloEngine.cpp.

178  {
179  SetNChains(i);
180  obval = new double[fMCMCNChains * kmax];
181  obweight = new double[fMCMCNChains * kwmax];
182 }

◆ setNoLegend()

void MonteCarloEngine::setNoLegend ( bool  legend)
inline

A set method to toggle the printing of legends in 1D and 2D histograms.

Parameters
[in]legenda boolean to toggle the printing of legends

Definition at line 307 of file MonteCarloEngine.h.

308  {
309  noLegend = legend;
310  };

◆ setPrintLoglikelihoodPlots()

void MonteCarloEngine::setPrintLoglikelihoodPlots ( bool  plot)
inline

A set method to toggle the printing of Parameters vs. Loglikelihood.

Parameters
[in]legenda boolean to toggle the printing of the 2D plots

Definition at line 316 of file MonteCarloEngine.h.

317  {
318  PrintLoglikelihoodPlots = plot;
319  };

◆ setPrintLogo()

void MonteCarloEngine::setPrintLogo ( bool  print)
inline

A set method to toggle the printing of logo on the histogram pdf.

Parameters
[in]printa boolean to toggle the printing

Definition at line 298 of file MonteCarloEngine.h.

299  {
300  printLogo = print;
301  };

◆ setSignificants()

void MonteCarloEngine::setSignificants ( unsigned int  significants_i)
inline

A set method to set the number of significant digits in the output.

Parameters
[in]significants_ithe number of significant digits in the output

Definition at line 404 of file MonteCarloEngine.h.

405  {
406  significants = significants_i;
407  };

◆ setSmooth()

void MonteCarloEngine::setSmooth ( int  int_N)
inline

A set method to set the number of smoothing passes for ROOT in 1D histograms.

Parameters
[in]int_Nnumber of smoothing passes for ROOT

Definition at line 368 of file MonteCarloEngine.h.

369  {
370  nSmooth = int_N;
371  };

◆ setWriteLogLikelihoodChain()

void MonteCarloEngine::setWriteLogLikelihoodChain ( bool  LL)
inline

A set method to toggle the writing of Loglikelihood in the ROOT tree.

Parameters
[in]LLa boolean to toggle the writing of Loglikelihood in the ROOT tree

Definition at line 325 of file MonteCarloEngine.h.

326  {
327  WriteLogLikelihoodChain = LL;
328  };

◆ setWriteParametersChain()

void MonteCarloEngine::setWriteParametersChain ( bool  LL)
inline

A set method to toggle the writing of parameters in the ROOT tree.

Parameters
[in]LLa boolean to toggle the writing of parameters in the ROOT tree

Definition at line 342 of file MonteCarloEngine.h.

343  {
344  WriteParametersChain = LL;
345  };

◆ writePreRunData()

std::string MonteCarloEngine::writePreRunData ( )

A method to write in a text file the best fit parameters and the prerun scale factors.

Returns
a string containing the data

Definition at line 1306 of file MonteCarloEngine.cpp.

1307 {
1308  std::vector<double> mode(GetBestFitParameters());
1309  if (mode.size() == 0) {
1310  throw std::runtime_error("\n ERROR: Global Mode could not be determined possibly because of infinite loglikelihood. PreRun Data cannot be stored.\n");
1311  }
1312  std::vector<double> scales(GetScaleFactors().at(0));
1313  std::ostringstream StatsLog;
1314  for (unsigned int i = 0; i < mode.size(); i++)
1315  StatsLog << GetParameter(i).GetName() << " " << mode.at(i) << " " << scales.at(i) << std::endl;
1316  return StatsLog.str().c_str();
1317 }

Member Data Documentation

◆ alpha2D

double MonteCarloEngine::alpha2D
private

A number between 0. and 1. that sets the opacity level of 2D Histograms, 1. being fully opaque.

Definition at line 484 of file MonteCarloEngine.h.

◆ CGO

std::vector<CorrelatedGaussianObservables>& MonteCarloEngine::CGO
private

A vector of correlated Gaussian observables.

Definition at line 445 of file MonteCarloEngine.h.

◆ CGP

const std::vector<CorrelatedGaussianParameters>& MonteCarloEngine::CGP
private

A vector of correlated Gaussian parameters.

Definition at line 442 of file MonteCarloEngine.h.

◆ cindex

unsigned int MonteCarloEngine::cindex
private

An index to distinguish between succesive canvases used to draw histograms.

Definition at line 474 of file MonteCarloEngine.h.

◆ CorrelationMap

std::map<std::string, TPrincipal *> MonteCarloEngine::CorrelationMap
private

A map between the name of a theory observable and its maximum value.

Definition at line 453 of file MonteCarloEngine.h.

◆ DPars

std::map<std::string, double> MonteCarloEngine::DPars
private

A map between parameter names and their values.

Definition at line 447 of file MonteCarloEngine.h.

◆ DPars_allChains

std::vector<std::map<std::string, double> > MonteCarloEngine::DPars_allChains
private

A vector of maps between parameter names and their values for all chains.

Definition at line 448 of file MonteCarloEngine.h.

◆ gIdx

int MonteCarloEngine::gIdx
private

Definition at line 485 of file MonteCarloEngine.h.

◆ HEPfit_green

TColor* MonteCarloEngine::HEPfit_green
private

< The number of significant digits in the Statistics File.

Definition at line 490 of file MonteCarloEngine.h.

◆ HEPfit_red

TColor* MonteCarloEngine::HEPfit_red
private

< The colour green for HEPfit.

Definition at line 491 of file MonteCarloEngine.h.

◆ Histo1D

std::map<std::string, BCH1D> MonteCarloEngine::Histo1D
private

A map between pointers to objects of type BCH1D (BAT) and their given names.

Definition at line 449 of file MonteCarloEngine.h.

◆ Histo2D

std::map<std::string, BCH2D> MonteCarloEngine::Histo2D
private

A map between pointers to objects of type BCH2D (BAT) and their given names.

Definition at line 450 of file MonteCarloEngine.h.

◆ histogram2Dtype

int MonteCarloEngine::histogram2Dtype
private

Type of 2D Histogram 1001 -> box pixel, 101 -> filled, 1 -> contour.

Definition at line 479 of file MonteCarloEngine.h.

◆ histogramBufferSize

unsigned int MonteCarloEngine::histogramBufferSize
private

< The colour red for HEPfit.

Definition at line 492 of file MonteCarloEngine.h.

◆ HistoLog

std::ostringstream MonteCarloEngine::HistoLog
private

A stream to store the output messages from printing and checking histograms.

Definition at line 459 of file MonteCarloEngine.h.

◆ hMCMCLogLikelihood

double MonteCarloEngine::hMCMCLogLikelihood
private

A variable containing the LogLikelihood values to be put into the ROOT tree.

Definition at line 471 of file MonteCarloEngine.h.

◆ hMCMCLogPriorProbability

double MonteCarloEngine::hMCMCLogPriorProbability
private

A variable containing the LogPriorProbability values to be put into the ROOT tree.

Definition at line 473 of file MonteCarloEngine.h.

◆ hMCMCLogProbability

double MonteCarloEngine::hMCMCLogProbability
private

A variable containing the LogProbability values to be put into the ROOT tree.

Definition at line 472 of file MonteCarloEngine.h.

◆ hMCMCObservables

std::vector<std::vector<double> > MonteCarloEngine::hMCMCObservables
private

A vector of vectors containing the observables values of all the chains to be put into the ROOT tree.

Definition at line 465 of file MonteCarloEngine.h.

◆ hMCMCObservables_weight

std::vector<std::vector<double> > MonteCarloEngine::hMCMCObservables_weight
private

A vector of vectors containing the observables weight of all the chains to be put into the ROOT tree.

Definition at line 466 of file MonteCarloEngine.h.

◆ hMCMCObservableTree

TTree* MonteCarloEngine::hMCMCObservableTree
private

A ROOT tree that contains the observables values and weight when the chains are written.

Definition at line 463 of file MonteCarloEngine.h.

◆ hMCMCParameters

std::vector<std::vector<double> > MonteCarloEngine::hMCMCParameters
private

A vector of vectors containing the parameter values of all the chains to be put into the ROOT tree.

Definition at line 467 of file MonteCarloEngine.h.

◆ hMCMCParameterTree

TTree* MonteCarloEngine::hMCMCParameterTree
private

A ROOT tree that contains the parameter values when the chains are written.

Definition at line 464 of file MonteCarloEngine.h.

◆ hMCMCTree_Observables

std::vector<double> MonteCarloEngine::hMCMCTree_Observables
private

A vector containing the observables values to be put into the ROOT tree.

Definition at line 468 of file MonteCarloEngine.h.

◆ hMCMCTree_Observables_weight

std::vector<double> MonteCarloEngine::hMCMCTree_Observables_weight
private

A vector containing the observables weight to be put into the ROOT tree.

Definition at line 469 of file MonteCarloEngine.h.

◆ hMCMCTree_Parameters

std::vector<double> MonteCarloEngine::hMCMCTree_Parameters
private

A vector containing the parameter values to be put into the ROOT tree.

Definition at line 470 of file MonteCarloEngine.h.

◆ kchainedObs

unsigned int MonteCarloEngine::kchainedObs
private

The number of observables for which the chains should be written.

Definition at line 458 of file MonteCarloEngine.h.

◆ kmax

unsigned int MonteCarloEngine::kmax
private

The number of observables.

Definition at line 457 of file MonteCarloEngine.h.

◆ kwmax

unsigned int MonteCarloEngine::kwmax
private

The number of observables whose weights are used for the MCMC.

Definition at line 456 of file MonteCarloEngine.h.

◆ LogLikelihood_max

double MonteCarloEngine::LogLikelihood_max
private

< vector to store the value of the parameters at maximum LogLikelihood;

Definition at line 494 of file MonteCarloEngine.h.

◆ Mod

StandardModel* MonteCarloEngine::Mod
private

A pointer to an abject of type Model.

Definition at line 446 of file MonteCarloEngine.h.

◆ ModPars

const std::vector<ModelParameter>& MonteCarloEngine::ModPars
private

A vector of model parameters.

Definition at line 441 of file MonteCarloEngine.h.

◆ nBins1D

unsigned int MonteCarloEngine::nBins1D
private

The number of bins in a 1D histogram.

Definition at line 487 of file MonteCarloEngine.h.

◆ nBins2D

unsigned int MonteCarloEngine::nBins2D
private

The number of bins in a 2D histogram.

Definition at line 488 of file MonteCarloEngine.h.

◆ noLegend

bool MonteCarloEngine::noLegend
private

A flag to toggle the histogram legends.

Definition at line 480 of file MonteCarloEngine.h.

◆ nSmooth

int MonteCarloEngine::nSmooth
private

The number of times a 1D histogram should be smoothed.

Definition at line 478 of file MonteCarloEngine.h.

◆ NumOfDiscardedEvents

int MonteCarloEngine::NumOfDiscardedEvents
private

The number of events for which the update of the model fails and these events are not used for the MCMC run.

Definition at line 461 of file MonteCarloEngine.h.

◆ NumOfUsedEvents

int MonteCarloEngine::NumOfUsedEvents
private

The number of events for which the model is successfully updated and hence used for the MCMC run.

Definition at line 460 of file MonteCarloEngine.h.

◆ Obs2D_ALL

std::vector<Observable2D>& MonteCarloEngine::Obs2D_ALL
private

A vector of all pairs of observable for Observable2D.

Definition at line 444 of file MonteCarloEngine.h.

◆ Obs_ALL

boost::ptr_vector<Observable>& MonteCarloEngine::Obs_ALL
private

A vector of all observables.

Definition at line 443 of file MonteCarloEngine.h.

◆ obval

double* MonteCarloEngine::obval
private

A pointer to the vector of observable values.

Definition at line 454 of file MonteCarloEngine.h.

◆ obweight

double* MonteCarloEngine::obweight
private

A pointer to the vector of observable weights.

Definition at line 455 of file MonteCarloEngine.h.

◆ ofi

std::ofstream MonteCarloEngine::ofi
private

Definition at line 475 of file MonteCarloEngine.h.

◆ par_at_LL_max

std::vector<double> MonteCarloEngine::par_at_LL_max
private

< The size of the buffer used for the histograms.

Definition at line 493 of file MonteCarloEngine.h.

◆ PrintLoglikelihoodPlots

bool MonteCarloEngine::PrintLoglikelihoodPlots
private

A flag to toggle the printing of Parameters vs. Loglikelihood.

Definition at line 481 of file MonteCarloEngine.h.

◆ printLogo

bool MonteCarloEngine::printLogo
private

A flag that is set to true for printing the logo on the histogram pdf.

Definition at line 477 of file MonteCarloEngine.h.

◆ rank

int MonteCarloEngine::rank
private

Rank of the process for a MPI run. Value is 0 for a serial run.

Definition at line 462 of file MonteCarloEngine.h.

◆ rIdx

int MonteCarloEngine::rIdx
private

Definition at line 486 of file MonteCarloEngine.h.

◆ significants

unsigned int MonteCarloEngine::significants
private

Definition at line 489 of file MonteCarloEngine.h.

◆ thMax

std::map<std::string, double> MonteCarloEngine::thMax
private

A map between the name of a theory observable and its maximum value.

Definition at line 452 of file MonteCarloEngine.h.

◆ thMin

std::map<std::string, double> MonteCarloEngine::thMin
private

A map between the name of a theory observable and its minimum value.

Definition at line 451 of file MonteCarloEngine.h.

◆ unknownParameters

std::vector<std::string> MonteCarloEngine::unknownParameters
private

A vector to contain the unkenown parameters passed in the configuration file.

Definition at line 476 of file MonteCarloEngine.h.

◆ WriteLogLikelihoodChain

bool MonteCarloEngine::WriteLogLikelihoodChain
private

A flag to toggle the writing of Loglikelihood chains in the ROOT tree.

Definition at line 482 of file MonteCarloEngine.h.

◆ WriteParametersChain

bool MonteCarloEngine::WriteParametersChain
private

A flag to toggle the writing of parameters chains in the ROOT tree.

Definition at line 483 of file MonteCarloEngine.h.

The documentation for this class was generated from the following files:
MonteCarloEngine::HEPfit_green
TColor * HEPfit_green
< The number of significant digits in the Statistics File.
Definition: MonteCarloEngine.h:490
Observable::computeTheoryValue
double computeTheoryValue()
A method to access the computed theory value of the observable.
Definition: Observable.cpp:130
MonteCarloEngine::hMCMCObservables
std::vector< std::vector< double > > hMCMCObservables
A vector of vectors containing the observables values of all the chains to be put into the ROOT tree.
Definition: MonteCarloEngine.h:465
MonteCarloEngine::ModPars
const std::vector< ModelParameter > & ModPars
A vector of model parameters.
Definition: MonteCarloEngine.h:441
MonteCarloEngine::WriteParametersChain
bool WriteParametersChain
A flag to toggle the writing of parameters chains in the ROOT tree.
Definition: MonteCarloEngine.h:483
MonteCarloEngine::nBins1D
unsigned int nBins1D
The number of bins in a 1D histogram.
Definition: MonteCarloEngine.h:487
MonteCarloEngine::HEPfit_red
TColor * HEPfit_red
< The colour green for HEPfit.
Definition: MonteCarloEngine.h:491
gslpp::matrix< double >
A class for constructing and defining operations on real matrices.
Definition: gslpp_matrix_double.h:48
MonteCarloEngine::PrintHistogram
void PrintHistogram(std::string &OutFile, Observable &it, const std::string OutputDir)
Overloaded from PrintHistogram(BCModelOutput&, const std::string) to print histogram for observables.
Definition: MonteCarloEngine.cpp:786
MonteCarloEngine::nBins2D
unsigned int nBins2D
The number of bins in a 2D histogram.
Definition: MonteCarloEngine.h:488
MonteCarloEngine::hMCMCTree_Observables_weight
std::vector< double > hMCMCTree_Observables_weight
A vector containing the observables weight to be put into the ROOT tree.
Definition: MonteCarloEngine.h:469
MonteCarloEngine::Obs_ALL
boost::ptr_vector< Observable > & Obs_ALL
A vector of all observables.
Definition: MonteCarloEngine.h:443
MonteCarloEngine::Obs2D_ALL
std::vector< Observable2D > & Obs2D_ALL
A vector of all pairs of observable for Observable2D.
Definition: MonteCarloEngine.h:444
MonteCarloEngine::alpha2D
double alpha2D
A number between 0. and 1. that sets the opacity level of 2D Histograms, 1. being fully opaque.
Definition: MonteCarloEngine.h:484
gslpp::log
complex log(const complex &z)
Definition: gslpp_complex.cpp:342
MonteCarloEngine::hMCMCObservableTree
TTree * hMCMCObservableTree
A ROOT tree that contains the observables values and weight when the chains are written.
Definition: MonteCarloEngine.h:463
MonteCarloEngine::Mod
StandardModel * Mod
A pointer to an abject of type Model.
Definition: MonteCarloEngine.h:446
MonteCarloEngine::hMCMCLogProbability
double hMCMCLogProbability
A variable containing the LogProbability values to be put into the ROOT tree.
Definition: MonteCarloEngine.h:472
MonteCarloEngine::Histo1D
std::map< std::string, BCH1D > Histo1D
A map between pointers to objects of type BCH1D (BAT) and their given names.
Definition: MonteCarloEngine.h:449
MonteCarloEngine::CheckHistogram
void CheckHistogram(TH1 &hist, const std::string name)
This member checks if there is overflow of the 1D histogram.
Definition: MonteCarloEngine.cpp:594
MonteCarloEngine::Print2D
void Print2D(BCH2D hist2D, const char *filename, int ww=0, int wh=0)
Definition: MonteCarloEngine.cpp:700
MonteCarloEngine::hMCMCTree_Parameters
std::vector< double > hMCMCTree_Parameters
A vector containing the parameter values to be put into the ROOT tree.
Definition: MonteCarloEngine.h:470
MonteCarloEngine::unknownParameters
std::vector< std::string > unknownParameters
A vector to contain the unkenown parameters passed in the configuration file.
Definition: MonteCarloEngine.h:476
MonteCarloEngine::thMin
std::map< std::string, double > thMin
A map between the name of a theory observable and its minimum value.
Definition: MonteCarloEngine.h:451
MonteCarloEngine::thMax
std::map< std::string, double > thMax
A map between the name of a theory observable and its maximum value.
Definition: MonteCarloEngine.h:452
MonteCarloEngine::InChainFillParametersTree
void InChainFillParametersTree()
Definition: MonteCarloEngine.cpp:964
MonteCarloEngine::getPrecision
int getPrecision(double value, double rms)
A get method to calculate the precision of the numbers in the Statistics file based on the rms.
Definition: MonteCarloEngine.cpp:1004
MonteCarloEngine::par_at_LL_max
std::vector< double > par_at_LL_max
< The size of the buffer used for the histograms.
Definition: MonteCarloEngine.h:493
MonteCarloEngine::histogramBufferSize
unsigned int histogramBufferSize
< The colour red for HEPfit.
Definition: MonteCarloEngine.h:492
MonteCarloEngine::getchainedObsSize
int getchainedObsSize() const
A get method to access the number of Observable chains requested.
Definition: MonteCarloEngine.h:195
MonteCarloEngine::hMCMCParameters
std::vector< std::vector< double > > hMCMCParameters
A vector of vectors containing the parameter values of all the chains to be put into the ROOT tree.
Definition: MonteCarloEngine.h:467
MonteCarloEngine::setDParsFromParameters
void setDParsFromParameters(const std::vector< double > &parameters, std::map< std::string, double > &DPars_i)
A method to rotate the diagonalized parameters to the original basis for correlated parameters.
Definition: MonteCarloEngine.cpp:261
MonteCarloEngine::significants
unsigned int significants
Definition: MonteCarloEngine.h:489
MonteCarloEngine::DPars
std::map< std::string, double > DPars
A map between parameter names and their values.
Definition: MonteCarloEngine.h:447
MonteCarloEngine::nSmooth
int nSmooth
The number of times a 1D histogram should be smoothed.
Definition: MonteCarloEngine.h:478
MonteCarloEngine::LogLikelihood
double LogLikelihood(const std::vector< double > &parameters)
This member calculates the loglikelihood for the observables included in the MCMC run.
Definition: MonteCarloEngine.cpp:312
gslpp::sqrt
complex sqrt(const complex &z)
Definition: gslpp_complex.cpp:385
MonteCarloEngine::DPars_allChains
std::vector< std::map< std::string, double > > DPars_allChains
A vector of maps between parameter names and their values for all chains.
Definition: MonteCarloEngine.h:448
MonteCarloEngine::LogLikelihood_max
double LogLikelihood_max
< vector to store the value of the parameters at maximum LogLikelihood;
Definition: MonteCarloEngine.h:494
MonteCarloEngine::gIdx
int gIdx
Definition: MonteCarloEngine.h:485
MonteCarloEngine::FirstDerivative
double FirstDerivative(BCParameter par, std::vector< double > point)
A method to calculate the first derivative.
Definition: MonteCarloEngine.cpp:1389
MonteCarloEngine::kchainedObs
unsigned int kchainedObs
The number of observables for which the chains should be written.
Definition: MonteCarloEngine.h:458
MonteCarloEngine::hMCMCLogLikelihood
double hMCMCLogLikelihood
A variable containing the LogLikelihood values to be put into the ROOT tree.
Definition: MonteCarloEngine.h:471
MonteCarloEngine::HistoLog
std::ostringstream HistoLog
A stream to store the output messages from printing and checking histograms.
Definition: MonteCarloEngine.h:459
MonteCarloEngine::rank
int rank
Rank of the process for a MPI run. Value is 0 for a serial run.
Definition: MonteCarloEngine.h:462
MonteCarloEngine::NumOfUsedEvents
int NumOfUsedEvents
The number of events for which the model is successfully updated and hence used for the MCMC run.
Definition: MonteCarloEngine.h:460
MonteCarloEngine::noLegend
bool noLegend
A flag to toggle the histogram legends.
Definition: MonteCarloEngine.h:480
MonteCarloEngine::rIdx
int rIdx
Definition: MonteCarloEngine.h:486
MonteCarloEngine::Histo2D
std::map< std::string, BCH2D > Histo2D
A map between pointers to objects of type BCH2D (BAT) and their given names.
Definition: MonteCarloEngine.h:450
MonteCarloEngine::DefineParameters
void DefineParameters()
A member to classify the prior of the model parameters varied in the Monte Carlo.
Definition: MonteCarloEngine.cpp:206
MonteCarloEngine::CorrelationMap
std::map< std::string, TPrincipal * > CorrelationMap
A map between the name of a theory observable and its maximum value.
Definition: MonteCarloEngine.h:453
StandardModel::Update
virtual bool Update(const std::map< std::string, double > &DPars)
The update method for StandardModel.
Definition: StandardModel.cpp:209
MonteCarloEngine::kmax
unsigned int kmax
The number of observables.
Definition: MonteCarloEngine.h:457
gslpp::log10
complex log10(const complex &z)
Definition: gslpp_complex.cpp:351
MonteCarloEngine::WriteLogLikelihoodChain
bool WriteLogLikelihoodChain
A flag to toggle the writing of Loglikelihood chains in the ROOT tree.
Definition: MonteCarloEngine.h:482
MonteCarloEngine::cindex
unsigned int cindex
An index to distinguish between succesive canvases used to draw histograms.
Definition: MonteCarloEngine.h:474
MonteCarloEngine::PrintLoglikelihoodPlots
bool PrintLoglikelihoodPlots
A flag to toggle the printing of Parameters vs. Loglikelihood.
Definition: MonteCarloEngine.h:481
MonteCarloEngine::hMCMCParameterTree
TTree * hMCMCParameterTree
A ROOT tree that contains the parameter values when the chains are written.
Definition: MonteCarloEngine.h:464
MonteCarloEngine::InChainFillObservablesTree
void InChainFillObservablesTree()
Definition: MonteCarloEngine.cpp:949
MonteCarloEngine::histogram2Dtype
int histogram2Dtype
Type of 2D Histogram 1001 -> box pixel, 101 -> filled, 1 -> contour.
Definition: MonteCarloEngine.h:479
MonteCarloEngine::CGO
std::vector< CorrelatedGaussianObservables > & CGO
A vector of correlated Gaussian observables.
Definition: MonteCarloEngine.h:445
Observable::getName
std::string getName() const
A get method to access the name of the observable.
Definition: Observable.h:323
MonteCarloEngine::kwmax
unsigned int kwmax
The number of observables whose weights are used for the MCMC.
Definition: MonteCarloEngine.h:456
MonteCarloEngine::NumOfDiscardedEvents
int NumOfDiscardedEvents
The number of events for which the update of the model fails and these events are not used for the MC...
Definition: MonteCarloEngine.h:461
MonteCarloEngine::CGP
const std::vector< CorrelatedGaussianParameters > & CGP
A vector of correlated Gaussian parameters.
Definition: MonteCarloEngine.h:442
MonteCarloEngine::SecondDerivative
double SecondDerivative(BCParameter par1, BCParameter par2, std::vector< double > point)
A method to calculate the second derivative.
Definition: MonteCarloEngine.cpp:1354
MonteCarloEngine::Print1D
void Print1D(BCH1D hist1D, const char *filename, int ww=0, int wh=0)
Definition: MonteCarloEngine.cpp:619
gslpp::exp
complex exp(const complex &z)
Definition: gslpp_complex.cpp:333
MonteCarloEngine::hMCMCObservables_weight
std::vector< std::vector< double > > hMCMCObservables_weight
A vector of vectors containing the observables weight of all the chains to be put into the ROOT tree.
Definition: MonteCarloEngine.h:466
MonteCarloEngine::obval
double * obval
A pointer to the vector of observable values.
Definition: MonteCarloEngine.h:454
MonteCarloEngine::Function_h
double Function_h(std::vector< double > point)
A method to calculate the LogLikelihood + LogAprioriProbability.
Definition: MonteCarloEngine.cpp:1424
QCD::getUnknownParameters
std::vector< std::string > getUnknownParameters()
A method to get the vector of the parameters that have been specified in the configuration file but n...
Definition: QCD.h:467
MonteCarloEngine::obweight
double * obweight
A pointer to the vector of observable weights.
Definition: MonteCarloEngine.h:455
MonteCarloEngine::hMCMCTree_Observables
std::vector< double > hMCMCTree_Observables
A vector containing the observables values to be put into the ROOT tree.
Definition: MonteCarloEngine.h:468
MonteCarloEngine::printLogo
bool printLogo
A flag that is set to true for printing the logo on the histogram pdf.
Definition: MonteCarloEngine.h:477
MonteCarloEngine::hMCMCLogPriorProbability
double hMCMCLogPriorProbability
A variable containing the LogPriorProbability values to be put into the ROOT tree.
Definition: MonteCarloEngine.h:473