BCModel.h

Go to the documentation of this file.

#ifndef __BCMODEL__H
#define __BCMODEL__H

// ---------------------------------------------------------

#include <vector>
#include <string>

#include "BCIntegrate.h"

// ROOT classes
class TNamed;
class TH1;
class TH2D;
class TGraph;
class TCanvas;
class TPostscript;
class TF1;

//BAT classes
class BCDataPoint;
class BCDataSet;
class BCParameter;
class BCH1D;
class BCH2D;

const int MAXNDATAPOINTVALUES = 20;

// ---------------------------------------------------------

class BCModel : public BCIntegrate
{

   public:

      BCModel();

      BCModel(const char * name);

      BCModel(const BCModel & bcmodel);

      virtual ~BCModel();

      BCModel & operator = (const BCModel & bcmodel);
      
      std::string GetName()
         { return fName; }

      int GetIndex()
         { return fIndex; }

      double GetModelAPrioriProbability()
         { return fModelAPriori; }

      double GetModelAPosterioriProbability()
         { return fModelAPosteriori; }

      double GetNormalization()
         { return fNormalization; }

      BCDataSet* GetDataSet()
         { return fDataSet; }

      BCDataPoint* GetDataPointLowerBoundaries()
         { return fDataPointLowerBoundaries; }

      BCDataPoint* GetDataPointUpperBoundaries()
         { return fDataPointUpperBoundaries; }

      double GetDataPointLowerBoundary(unsigned int index)
         { return fDataPointLowerBoundaries -> GetValue(index); }

      double GetDataPointUpperBoundary(unsigned int index)
         { return fDataPointUpperBoundaries -> GetValue(index); }

      bool GetFlagBoundaries();

      int GetNDataPoints();

      BCDataPoint * GetDataPoint(unsigned int index);

      unsigned int GetNDataPointsMinimum()
         { return fNDataPointsMinimum; }

      unsigned int GetNDataPointsMaximum()
         { return fNDataPointsMaximum; }

      unsigned int GetNParameters()
         { return fParameterSet ? fParameterSet -> size() : 0; }

      BCParameter * GetParameter(int index);

      BCParameter * GetParameter(const char * name);

      BCParameterSet * GetParameterSet()
         { return fParameterSet; }

      double GetBestFitParameter(unsigned int index);

      double GetBestFitParameterError(unsigned int index);

      std::vector <double> GetBestFitParameters()
         { return fBestFitParameters; }
      
      std::vector <double> GetBestFitParameterErrors() 
         { return fBestFitParameterErrors; }

      double GetBestFitParameterMarginalized(unsigned int index);

      std::vector <double> GetBestFitParametersMarginalized()
         { return fBestFitParametersMarginalized; }

      TH2D * GetErrorBandXY()
         { return fErrorBandXY; }

      TH2D * GetErrorBandXY_yellow(double level=.68, int nsmooth=0);

      std::vector <double> GetErrorBand(double level);

      TGraph * GetErrorBandGraph(double level1, double level2);

      TGraph * GetFitFunctionGraph(std::vector <double> parameters);

      TGraph * GetFitFunctionGraph()
         { return this -> GetFitFunctionGraph(this -> GetBestFitParameters()); }

      TGraph * GetFitFunctionGraph(std::vector <double> parameters, double xmin, double xmax, int n=1000);

      bool GetFixedDataAxis(unsigned int index);

      void SetName(const char * name)
         { fName = name; }

      void SetIndex(int index)
         { fIndex = index; }

      void SetParameterSet( BCParameterSet * parset )
         { fParameterSet = parset; }

      int SetParameterRange(int index, double parmin, double parmax);

      void SetModelAPrioriProbability(double probability)
         { fModelAPriori = probability; }

      void SetModelAPosterioriProbability(double probability)
         { fModelAPosteriori = probability; }

      void SetNormalization(double norm)
         { fNormalization = norm; }

      void SetDataSet(BCDataSet* dataset)
         { fDataSet = dataset; fNormalization = -1.0; }

      void SetSingleDataPoint(BCDataPoint * datapoint);

      void SetSingleDataPoint(BCDataSet * dataset, unsigned int index);

      void SetNDataPointsMinimum(unsigned int minimum)
         { fNDataPointsMinimum = minimum; }

      void SetNDataPointsMaximum(unsigned int maximum)
         { fNDataPointsMaximum = maximum; }

      void SetDataBoundaries(unsigned int index, double lowerboundary, double upperboundary, bool fixed=false);

      void SetErrorBandContinuous(bool flag);

      void SetNbins(const char * parname, int nbins);

      int SetPrior(int index, TF1* f); 

      int SetPrior(const char* name, TF1* f); 

      int SetPriorDelta(int index, double value); 

      int SetPriorDelta(const char* name, double value); 

      int SetPriorGauss(int index, double mean, double sigma); 

      int SetPriorGauss(const char* name, double mean, double sigma); 

      int SetPriorGauss(int index, double mean, double sigmadown, double sigmaup);

      int SetPriorGauss(const char* name, double mean, double sigmadown, double sigmaup);

      int SetPrior(int index, TH1 * h, bool flag=false); 

      int SetPrior(const char* name, TH1 * h, bool flag=false); 

      int SetPriorConstant(int index);

      int SetPriorConstant(const char* name);

      int SetPriorConstantAll();

      int AddParameter(const char * name, double lowerlimit, double upperlimit);

      int AddParameter(BCParameter* parameter);

      double APrioriProbability(std::vector <double> parameters)
         { return exp( this->LogAPrioriProbability(parameters) ); }

      virtual double LogAPrioriProbability(std::vector <double> parameters); 

      double Likelihood(std::vector <double> parameter)
         { return exp( this->LogLikelihood(parameter) ); }

      virtual double LogLikelihood(std::vector <double> parameter);

      double ProbabilityNN(std::vector <double> parameter)
         { return exp( this->LogProbabilityNN(parameter) ); }

      double LogProbabilityNN(std::vector <double> parameter);

      double Probability(std::vector <double> parameter)
         { return exp( this->LogProbability(parameter) ); }

      double LogProbability(std::vector <double> parameter);

      double ConditionalProbabilityEntry(BCDataPoint * datapoint, std::vector <double> parameters)
         { return exp( this->LogConditionalProbabilityEntry(datapoint, parameters) ); }

      virtual double LogConditionalProbabilityEntry(BCDataPoint * /*datapoint*/, std::vector <double> /*parameters*/)
         { flag_ConditionalProbabilityEntry = false; return 0.; }

      virtual double SamplingFunction(std::vector <double> parameters);

      double Eval(std::vector <double> parameters)
         { return exp( this->LogEval(parameters) ); }

      double LogEval(std::vector <double> parameters);

      double EvalSampling(std::vector <double> parameters);

      double Normalize();

      int CheckParameters(std::vector <double> parameters);

      void FindMode(std::vector<double> start = std::vector<double>(0));

      void FindModeMinuit(std::vector<double> start = std::vector<double>(0), int printlevel = 1);

      void WriteMode(const char * file);

      int ReadMode(const char * file);

      int ReadMarginalizedFromFile(const char * file);

      int ReadErrorBandFromFile(const char * file);

      int MarginalizeAll();

      BCH1D * GetMarginalized(BCParameter * parameter);

      BCH1D * GetMarginalized(const char * name)
         { return this -> GetMarginalized(this -> GetParameter(name)); }

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

      BCH2D * GetMarginalized(const char * name1, const char * name2)
         { return this -> GetMarginalized(this -> GetParameter(name1), this -> GetParameter(name2)); }

      int PrintAllMarginalized1D(const char * filebase);
      int PrintAllMarginalized2D(const char * filebase);
      int PrintAllMarginalized(const char * file, unsigned int hdiv=1, unsigned int ndiv=1);

      virtual void CorrelateDataPointValues(std::vector<double> &x);

      double GetPvalueFromChi2(std::vector<double> par, int sigma_index);

      double GetPvalueFromChi2Johnson(std::vector<double> par);

      double GetPvalueFromKolmogorov(const std::vector<double>& par, int index);


      double GetChi2Johnson(std::vector<double> par, const int nBins=-1);

      double GetAvalueFromChi2Johnson(TTree* tree, TH1D* distribution=0);

      double GetPvalueFromChi2NDoF(std::vector<double> par, int sigma_index);

      BCH1D * CalculatePValue(std::vector<double> par, bool flag_histogram=false);

      double GetPValue()
         { return fPValue; }

      double GetPValueNDoF()
         { return fPValueNDoF; }

      double GetChi2NDoF()
         { return fChi2NDoF; }

      std::vector<double> GetChi2Runs(int dataIndex, int sigmaIndex);

      void SetGoFNIterationsMax(int n)
         { fGoFNIterationsMax=n; }

      void SetGoFNIterationsRun(int n)
         { fGoFNIterationsRun=n; }

      void SetGoFNChains(int n)
         { fGoFNChains=n; }

      double HessianMatrixElement(BCParameter * parameter1, BCParameter * parameter2, std::vector<double> point);

      void PrintSummary();

      void PrintResults(const char * file);

      void PrintShortFitSummary(int chi2flag=0);

      void PrintHessianMatrix(std::vector<double> parameters);

      void FixDataAxis(unsigned int index, bool fixed);


      virtual double CDF(const std::vector<double>& /*parameters*/,  int /*index*/, bool /*lower=false*/)
      {return 0.0;}


      int ResetResults();

   protected:

      int fIndex;

      std::string fName;

      double fModelAPriori;

      double fModelAPosteriori;

      BCParameterSet * fParameterSet;

      BCDataSet * fDataSet;

      unsigned int fNDataPointsMinimum;

      unsigned int fNDataPointsMaximum;

      bool flag_ConditionalProbabilityEntry;

      double fPValue;

      double fChi2NDoF;
      double fPValueNDoF;

      bool flag_discrete;

      int fGoFNIterationsMax;

      int fGoFNIterationsRun;

      int fGoFNChains;

//    std::vector<TF1*> fPriorContainer;
      std::vector<TNamed*> fPriorContainer;

      bool fPriorConstantAll;

      double fPriorConstantValue;

      std::vector<bool> fPriorContainerConstant;

      std::vector<bool> fPriorContainerInterpolate;

   private:

      double fNormalization;

      int CompareStrings(const char * string1, const char * string2);

      int NumberBins()
         { return (int) exp( 0.4 * log( (float)GetNDataPoints() ) ) + 2; }

      void RecalculatePriorConstant();

      void StoreMode();

      BCDataPoint * VectorToDataPoint(std::vector<double> data);
};

// ---------------------------------------------------------

typedef std::vector<BCModel*> BCModelContainer;

// ---------------------------------------------------------

#endif