MonteCarloEngine.h

Go to the documentation of this file.

/*
 * Copyright (C) 2012 HEPfit Collaboration
 * All rights reserved.
 *
 * 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 <BAT/BCModelOutput.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TFile.h>
#include <TRandom3.h>
#include <TPrincipal.h>
#include <map>
#include <string>
#include <sstream>
#include <boost/ptr_container/ptr_vector.hpp>

#define NBINSMODELPARS 100
#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(Model* Mod_i);

    void DefineParameters();

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

    void CheckHistogram(const TH1D& hist, const std::string name);

    void CheckHistogram(const TH2D& hist, const std::string name);

    void PrintHistogram(BCModelOutput & out, Observable & it, const std::string OutputDir);

    void PrintHistogram(BCModelOutput& out, const std::string OutputDir);

    void MCMCIterationInterface();

    void setNChains(unsigned int i);

    void AddChains();

    std::string getHistoLog() const
    {
        return HistoLog.str().c_str();
    }

    std::string computeStatistics();

    std::string writePreRunData();

    void PrintCorrelationMatrix(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;
    }
    
    void setHistogramOverFlow(bool overflow)
    {
        if (overflow) TH1::StatOverflows(kTRUE);
        else TH1::StatOverflows(kFALSE);
    }
    
    double computeNormalization();
    
    double SecondDerivative(const BCParameter * par1, const BCParameter * par2, std::vector<double> point);
    
    double FirstDerivative(const 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);
    
private:
    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; 
    Model* Mod; 
    std::map<std::string, double> DPars; 
    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; 
    std::ostringstream HistoLog; 
    int NumOfUsedEvents; 
    int NumOfDiscardedEvents; 
    int rank; 

};

#endif