MonteCarloEngine.h

Go to the documentation of this file.

/*
 * Copyright (C) 2012 HEPfit Collaboration
 *
 *
 * For the licensing terms see doc/COPYING.
 */
 
#ifndef __MONTECARLOENGINE__H
#define __MONTECARLOENGINE__H
 
#include "Observable.h"
#include "Observable2D.h"
#include "CorrelatedGaussianObservables.h"
#include "CorrelatedGaussianParameters.h"
#include "ModelParameter.h"
#include "Model.h"
#include <BAT/BCModel.h>
#include <BAT/BCH1D.h>
#include <BAT/BCH2D.h>
#include <TFile.h>
#include <TPrincipal.h>
#include <TColor.h>
#include <map>
#include <string>
#include <sstream>
#include <boost/ptr_container/ptr_vector.hpp>
 
#define NBINS1D 100
#define NBINS2D 100
#define NSTEPS 1e5
 
class MonteCarloEngine : public BCModel {
public:
 
    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);
 
    ~MonteCarloEngine();
 
    void Initialize(StandardModel* Mod_i);
    
    void CreateHistogramMaps();
 
    void DefineParameters();
 
    double LogLikelihood(const std::vector <double>& parameters);
 
    void CheckHistogram(TH1& hist, const std::string name);
 
    void CheckHistogram(TH2& hist, const std::string name);
    
    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 PrintHistogram(std::string& OutFile, Observable & it, const std::string OutputDir);
 
    void PrintHistogram(std::string& OutFile, const std::string OutputDir);
 
    void MCMCUserIterationInterface();
 
    void setNChains(unsigned int i);
 
    void AddChains();
    
    void InChainFillObservablesTree();
    
    void InChainFillParametersTree();
 
    std::string getHistoLog() const
    {
        return HistoLog.str().c_str();
    }
    
    int getchainedObsSize() const
    {
        return kchainedObs;
    }
 
    std::string computeStatistics();
 
    std::string writePreRunData();
 
    void PrintCorrelationMatrixToLaTeX(const std::string filename);
 
    int getNumOfDiscardedEvents() const
    {
        return NumOfDiscardedEvents;
    }
 
    int getNumOfUsedEvents() const
    {
        return NumOfUsedEvents;
    }
    
    std::map<std::string, BCH1D> getHistograms1D() const
    {
        return Histo1D;
    }
 
    std::map<std::string, BCH2D> getHistograms2D() const
    {
        return Histo2D;
    }
    
    double computeNormalizationLME();
    
    std::vector<double> computeNormalizationMC(int NIterationNormalizationMC);
    
    double SecondDerivative(BCParameter par1, BCParameter par2, std::vector<double> point);
    
    double FirstDerivative(BCParameter par, std::vector<double> point);
    
    double Function_h(std::vector<double> point);
    
    void setDParsFromParameters(const std::vector<double>& parameters, std::map<std::string,double>& DPars_i);
    
    void setPrintLogo(bool print)
    {
        printLogo = print;
    };
    
    void setNoLegend(bool legend)
    {
        noLegend = legend;
    };
    
    void setPrintLoglikelihoodPlots(bool plot)
    {
        PrintLoglikelihoodPlots = plot;
    };
    
    void setWriteLogLikelihoodChain(bool LL)
    {
        WriteLogLikelihoodChain = LL;
    };
    
    bool getWriteLogLikelihoodChain()
    {
        return WriteLogLikelihoodChain;
    };
    
    void setWriteParametersChain(bool LL)
    {
        WriteParametersChain = LL;
    };
    
    bool getWriteParametersChain()
    {
        return WriteParametersChain;
    };
    
    void setAlpha2D(double alpha)
    {
        alpha2D = alpha;
    };
    
    void setSmooth(int int_N)
    {
        nSmooth = int_N;
    };
    
    void setHistogram2DType(int type)
    {
        histogram2Dtype = type;
    };
    
    void setNBins1D(unsigned int nbins)
    {
        nBins1D = nbins;
    };
    
    void setNBins2D(unsigned int nbins)
    {
        nBins2D = nbins;
    };
    
    void setSignificants(unsigned int significants_i)
    {
        significants = significants_i;
    };
    
    void setHistogramBufferSize(unsigned int histogramBufferSize_i)
    {
        histogramBufferSize = histogramBufferSize_i;
    };
    
    int getMPIWorldSize() 
    {
#ifdef _MPI
        return procnum;
#else
        return 1;
#endif
    }
    
private:
    
    int getPrecision(double value, double rms);
    
    const std::vector<ModelParameter>& ModPars; 
    const std::vector<CorrelatedGaussianParameters>& CGP; 
    boost::ptr_vector<Observable>& Obs_ALL; 
    std::vector<Observable2D>& Obs2D_ALL; 
    std::vector<CorrelatedGaussianObservables>& CGO; 
    StandardModel* Mod; 
    std::map<std::string, double> DPars; 
    std::vector<std::map<std::string, double> > DPars_allChains; 
    std::map<std::string, BCH1D> Histo1D; 
    std::map<std::string, BCH2D> Histo2D; 
    std::map<std::string, double> thMin; 
    std::map<std::string, double> thMax; 
    std::map<std::string, TPrincipal *> CorrelationMap; 
    double *obval; 
    double *obweight; 
    unsigned int kwmax; 
    unsigned int kmax; 
    unsigned int kchainedObs; 
    std::ostringstream HistoLog; 
    int NumOfUsedEvents; 
    int NumOfDiscardedEvents; 
    int rank; 
    TTree * hMCMCObservableTree; 
    TTree * hMCMCParameterTree; 
    std::vector<std::vector<double> > hMCMCObservables; 
    std::vector<std::vector<double> > hMCMCObservables_weight; 
    std::vector<std::vector<double> > hMCMCParameters; 
    std::vector<double> hMCMCTree_Observables; 
    std::vector<double> hMCMCTree_Observables_weight; 
    std::vector<double> hMCMCTree_Parameters; 
    double hMCMCLogLikelihood; 
    double hMCMCLogProbability; 
    double hMCMCLogPriorProbability; 
    unsigned int cindex;
    std::ofstream ofi;
    std::vector<std::string> unknownParameters; 
    bool printLogo; 
    int nSmooth; 
    int histogram2Dtype; 
    bool noLegend; 
    bool PrintLoglikelihoodPlots; 
    bool WriteLogLikelihoodChain; 
    bool WriteParametersChain; 
    double alpha2D;
    int gIdx;
    int rIdx;
    unsigned int nBins1D; 
    unsigned int nBins2D; 
    unsigned int significants; 
    TColor * HEPfit_green; 
    TColor * HEPfit_red; 
    unsigned int histogramBufferSize; 
    std::vector<double> par_at_LL_max; 
    double LogLikelihood_max; 
};
 
#endif