BCTemplateFitter.h

Go to the documentation of this file.

#ifndef __BCTEMPLATEFITTER__H
#define __BCTEMPLATEFITTER__H

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

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

#include <vector>

#include <TH1D.h>
#include <TF1.h>

#include "../../BAT/BCModel.h"

class TH2D;

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

class BCTemplateFitter : public BCModel
{
   public:

      BCTemplateFitter();

      BCTemplateFitter(const char * name);

      ~BCTemplateFitter();

      int GetNTemplates()
         { return int(fTemplateParIndexContainer.size()); }

      int GetNTemplatesType(int templatetype);

      int GetNSystErrors()
         { return int(fSystErrorParIndexContainer.size()); }

      int GetNDF()
         { return fHistData.GetNbinsX() - GetNParameters(); }

      int GetNRatios()
         { return int(fHistRatios1D.size()); }

      std::vector<TH1D> GetHistRatios1D()
         { return fHistRatios1D; }

      std::vector<TH1D> GetTemplateHistogramContainer()
         { return fTemplateHistogramContainer; };

      std::vector<TF1> GetTemplateFunctionContainer()
         { return fTemplateFunctionContainer; };

      TH1D GetHistRatio1D(int index)
         { return fHistRatios1D.at(index); }

      TH2D* GetUncertaintyHistogramExp()
      { return fUncertaintyHistogramExp; };

      TH2D* GetUncertaintyHistogramObsPosterior()
      { return fUncertaintyHistogramObsPosterior; };

      int GetIndexTemplate(const char * name);

      int GetIndexSystError(const char * name);

      int GetParIndexTemplate(const char * name);

      int GetParIndexTemplate(int index);

      int GetParIndexEff(const char * name);

      int GetParIndexSystError(const char * name);

      //      TH1D GetTemplate(const char * name);

      TH1D GetData()
         { return fHistData; }

      int GetHistogramIndex(int templatenumber)
         { return fHistogramIndexContainer.at(templatenumber); };

      int GetFunctionIndex(int templatenumber)
         { return fFunctionIndexContainer.at(templatenumber); };

      std::vector<int> GetTemplateParIndexContainer()
         { return fTemplateParIndexContainer; };

      double Expectation(int binnumber, const std::vector<double>& parameters);

      double TemplateEfficiency(int templatenumber, int binnumber, const std::vector<double>& parameters);

      double TemplateProbability(int templatenumber, int binnumber, const std::vector<double>& parameters);


      void SetFlagPhysicalLimits(bool flag)
         { fFlagPhysicalLimits = flag; }

      void SetFlagFixNorm(bool flag)
         { fFlagFixNorm = flag; }

      void SetNorm(double norm)
         { fNorm = norm; }

      int SetData(const TH1D& hist);

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

      int Initialize();

      int AddTemplate(TH1D hist, const char * name, double Nmin=0, double Nmax=0);

      int AddTemplate(TF1 func, const char * name, double Nmin=0, double Nmax=0, const char* options="");

      int SetTemplateEfficiency(const char * name, double effmean = 1., double effsigma = 0.);

      int SetTemplateEfficiency(const char * name, TH1D eff, TH1D efferr);

      int AddSystError(const char * errorname, const char * errtype = "gauss");

      int SetTemplateSystError(const char * errorname, const char * templatename, TH1D parerror);

      //      int AddSystErrorCorrelation(const char * errorname1, const char * errnorame2, double corr);

      int ConstrainSum(std::vector<int> indices, double mean, double rms);

      int CalculateRatio(int index, std::vector<int> indices, double rmin = -1.0, double rmax = 1.0);

      int CompareHistogramProperties(TH1D hist1, TH1D hist2);

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

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

      double CalculateMaxLike();

      double CalculatePValue();

      double CalculateKSProb();

      void MCMCUserIterationInterface();

      void PrintStack(const char * filename = "stack.ps", const char * options="LE2E0D");

      void PrintRatios(const char * filename = "ratio.ps", int option = 0, double ovalue = 0.);

      int PrintTemplate(const char * name, const char * filename);

      void PrintTemp();

      TH1D CreateErrorHist(TH1D hist, TH1D histerr);

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

      int PerformFit();

      TH1D CombineUncertainties(const char * name);

    protected:

      TH1D fHistData;

      // histogram container

      std::vector<TH1D> fTemplateHistogramContainer;

      std::vector<TF1> fTemplateFunctionContainer;

      std::vector<TH1D> fEffHistogramContainer;

      std::vector<TH1D> fEffErrHistogramContainer;

      std::vector< std::vector<TH1D> > fSystErrorHistogramContainer;

      // name container

      std::vector<std::string> fTemplateNameContainer;

      std::vector<std::string> fSystErrorNameContainer;

      // index container

      std::vector<int> fTemplateParIndexContainer;

      std::vector<int> fHistogramIndexContainer;

      std::vector<int> fFunctionIndexContainer;

      std::vector<int> fEffParIndexContainer;

      std::vector<int> fSystErrorParIndexContainer;

      // error type container

      std::vector<std::string> fSystErrorTypeContainer;


      std::vector<int> fTemplateTypeContainer;

      std::vector< std::vector<int> > fConstraintSumIndices;

      std::vector< double > fConstraintSumMean;

      std::vector< double > fConstraintSumRMS;

      TH1D fHistNorm;

      std::vector< std::vector<int> > fIndicesRatios1D;

      std::vector<TH1D> fHistRatios1D;

      bool fFlagFixNorm;

      bool fFlagPhysicalLimits;

      TH2D * fUncertaintyHistogramExp;

      TH2D * fUncertaintyHistogramObsPosterior;

      double fNorm;

      int fNBins;

      double fXmin;

      double fXmax;

};

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

#endif