BCModel.h

Go to the documentation of this file.

#ifndef __BCMODEL__H
#define __BCMODEL__H

/*
 * Copyright (C) 2008, Daniel Kollar and Kevin Kroeninger.
 * All rights reserved.
 *
 * For the licensing terms see doc/COPYING.
 */

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

#include <vector>
#include <string>

#include "BAT/BCIntegrate.h"

// ROOT classes
class TH2D;
class TGraph;
class TCanvas;
class TPostscript;

//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);

      virtual ~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); };

      /*
       * Checks if the boundaries have been defined
       * @return true, if the boundaries have been set, false otherwise */
      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)
         { return fBestFitParameters[index]; };

      double GetBestFitParameterError(unsigned int index)
      { return fBestFitParameterErrors[index]; }; 

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

      double GetBestFitParameterMarginalized(unsigned int index)
         { return fBestFitParametersMarginalized[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; };

      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 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)
         { return 0.; };

      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 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);

      /*
       * @return The p-value */
      double GetPValue()
         { return fPValue; };

      double GetPValueNDoF()
         { return fPValueNDoF; };

      double GetChi2NDoF()
         { return fChi2NDoF; };

      /*
       * Set maximum number of iterations in the MCMC pre-run of the p-value
       * evaluation using MCMC */
      void SetGoFNIterationsMax(int n)
         { fGoFNIterationsMax=n; };

      /*
       * Set number of iterations in the MCMC normal run of the p-value
       * evaluation using MCMC */
      void SetGoFNIterationsRun(int n)
         { fGoFNIterationsRun=n; };

      /*
       * Set number of chains in the MCMC of the p-value
       * evaluation using MCMC */
      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;}


   /* @} */

   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;

      /*
       * Maximum number of iterations in the MCMC pre-run of the p-value
       * evaluation using MCMC */
      int fGoFNIterationsMax;

      /*
       * Number of iterations in the MCMC normal run of the p-value
       * evaluation using MCMC */
      int fGoFNIterationsRun;

      /*
       * Number of chains in the MCMC of the p-value
       * evaluation using MCMC */
      int fGoFNChains;

   private:

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

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

      double fNormalization;

      int NumberBins(){
         return (int)(exp(0.4 * log(this -> GetNDataPoints())) + 2);
      }

};

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

typedef std::vector<BCModel*> BCModelContainer;

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

#endif

---