BCTemplateFitter Class Reference

DEPRECETED class for fitting several templates to a data set. More...

#include <BCTemplateFitter.h>

Inherits BCModel.

Collaboration diagram for BCTemplateFitter:

List of all members.

Public Member Functions
	BCTemplateFitter ()
	BCTemplateFitter (const char *name)
	~BCTemplateFitter ()
int	GetNTemplates ()
int	GetNTemplatesType (int templatetype)
int	GetNSystErrors ()
int	GetNDF ()
int	GetNRatios ()
std::vector< TH1D >	GetHistRatios1D ()
std::vector< TH1D >	GetTemplateHistogramContainer ()
std::vector< TF1 >	GetTemplateFunctionContainer ()
TH1D	GetHistRatio1D (int index)
TH2D *	GetUncertaintyHistogramExp ()
TH2D *	GetUncertaintyHistogramObsPosterior ()
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	GetData ()
int	GetHistogramIndex (int templatenumber)
int	GetFunctionIndex (int templatenumber)
std::vector< int >	GetTemplateParIndexContainer ()
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)
void	SetFlagFixNorm (bool flag)
void	SetNorm (double 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	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 Attributes
TH1D	fHistData
std::vector< TH1D >	fTemplateHistogramContainer
std::vector< TF1 >	fTemplateFunctionContainer
std::vector< TH1D >	fEffHistogramContainer
std::vector< TH1D >	fEffErrHistogramContainer
std::vector< std::vector< TH1D > >	fSystErrorHistogramContainer
std::vector< std::string >	fTemplateNameContainer
std::vector< std::string >	fSystErrorNameContainer
std::vector< int >	fTemplateParIndexContainer
std::vector< int >	fHistogramIndexContainer
std::vector< int >	fFunctionIndexContainer
std::vector< int >	fEffParIndexContainer
std::vector< int >	fSystErrorParIndexContainer
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

---

Detailed Description

DEPRECETED class for fitting several templates to a data set.

This class is DEPRECETED and it will be removed in the future version of BAT. Use BCMTF instead.

This class can be used for fitting several template histograms to a data histogram. The templates are assumed to have no statistical uncertainty whereas the data are assumed to have Poissonian fluctuations in each bin. Several methods to judge the validity of the model are available.

Author:

Daniel Kollar

Kevin Kröninger

Date:

10.04.2010

Definition at line 40 of file BCTemplateFitter.h.

---

Constructor & Destructor Documentation

BCTemplateFitter::BCTemplateFitter

(

)

The default constructor.

Definition at line 32 of file BCTemplateFitter.cxx.

   : BCModel()
   , fFlagFixNorm(false)
   , fFlagPhysicalLimits(true)
   , fUncertaintyHistogramExp(0)
   , fUncertaintyHistogramObsPosterior(0)
   , fNorm(-1)
   , fNBins(-1)
   , fXmin(1.)
   , fXmax(0.)
{
   BCLog::OutWarning("Class BCTemplateFitter is depreceted and it will be removed");
   BCLog::OutWarning("in the future version of BAT. Use the new class BCMTF instead.");
}

BCTemplateFitter::BCTemplateFitter

(

const char *

name

)

A constructor.

Definition at line 48 of file BCTemplateFitter.cxx.

   : BCModel(name)
   , fFlagFixNorm(false)
   , fFlagPhysicalLimits(true)
   , fUncertaintyHistogramExp(0)
   , fUncertaintyHistogramObsPosterior(0)
   , fNorm(-1)
   , fNBins(-1)
   , fXmin(1.)
   , fXmax(0.)
{
   BCLog::OutWarning("Class BCTemplateFitter is depreceted and it will be removed");
   BCLog::OutWarning("in the future version of BAT. Use the new class BCMTF instead.");
}

BCTemplateFitter::~BCTemplateFitter

(

)

The default destructor.

Definition at line 64 of file BCTemplateFitter.cxx.

{
   if (fUncertaintyHistogramExp)
      delete fUncertaintyHistogramExp;

   if (fUncertaintyHistogramObsPosterior)
      delete fUncertaintyHistogramObsPosterior;
}

---

Member Function Documentation

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

Add a source of systematic uncertainty.

Parameters:

	errorname	The name of the source.
	errortype	The shape of the uncertainty in each bin.

Returns:

An error code.

Definition at line 406 of file BCTemplateFitter.cxx.

{
   // define parameter range
   double dx = 1.0;

   // check error type
   if (std::string(errtype) == std::string("gauss"))
      dx = 5.0;
   else if (std::string(errtype) == std::string("flat"))
      dx = 1.0;
   else {
      BCLog::OutError("BCTemplateFitter::AddSystError : Unknown error type.");
      return 0;
   }

   // add a parameter for the expectation value
   AddParameter(Form("%s", errorname), -dx, dx);

   // set prior for systematics
   SetPriorGauss(Form("%s", errorname), 0., 1.);

   // add name and index to containers
   fSystErrorNameContainer.push_back(errorname);
   fSystErrorParIndexContainer.push_back(GetNParameters()-1);

   // create histogram
   TH1D hist = TH1D(fHistData);

   // fill histograms
   for (int i = 1; i <= fNBins; ++i)
      hist.SetBinContent(i, 0.);

   // get number of templates
   int n = GetNTemplates();

   // define vector of histograms
   std::vector<TH1D> histvector;

   // add histograms
   for (int i = 0; i < n; ++i)
      histvector.push_back(hist);

   // add histogram vector
   fSystErrorHistogramContainer.push_back(histvector);

   // add error type to container
   fSystErrorTypeContainer.push_back(std::string(errtype));

   // no error
   return 1;
}

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

Add a template histogram. The templates do not have to be normalized. The histogram has to have the same number of bins and cover the same region as the data histogram.

Parameters:

	hist	The template histogram
	name	The name of the template
	Nmin	The lower limit of the normalization.
	Nmax	The upper limit of the normalization.

Definition at line 257 of file BCTemplateFitter.cxx.

{
   // check if histogram if filled
   if (hist.Integral() <= 0.) {
      BCLog::OutError("BCTemplateFitter::AddTemplate : Normalization is zero or less than that.");
      return 0;
   }

   // compare template properties with data
   if (fHistData.Integral() > 0.) {
      if (CompareHistogramProperties(fHistData, hist) != 1) {
         BCLog::OutError("BCTemplateFitter::AddTemplate : Data and template histogram properties are incompatible.");
         return 0;
      }
   }
   else {
      SetData( TH1D("", "", hist.GetNbinsX(), hist.GetXaxis()->GetXmin(), hist.GetXaxis()->GetXmax()) );
   }

   // check if prior makes sense
   if (fFlagPhysicalLimits && Nmin < 0)
      Nmin = 0;

   if (Nmin > Nmax) {
      BCLog::OutError("BCTemplateFitter::AddTemplate : Lower limit exceeds upper limit.");
      return 0;
   }

   // get number of templates
   int ntemplates = GetNTemplates();

   // set histogram color and style
   hist.SetFillColor(2 + ntemplates);
   hist.SetFillStyle(1001);
   hist.SetStats(kFALSE);

   // scale histogram
   hist.Scale(1.0 / hist.Integral());

   // check if template is consistent with other templates
   if (ntemplates > 0)
      if (!CompareHistogramProperties(fTemplateHistogramContainer.at(0), hist)) {
         BCLog::OutError("BCTemplateFitter::AddTemplate : Properties of template histogram is not compatible with older template histograms.");
         return 0;
      }

   // histograms
   TH1D histsysterror = TH1D(fHistData);

   // fill histograms
   for (int i = 1; i <= fNBins; ++i)
      histsysterror.SetBinContent(i, 0.);

   // get parameter index
   int parindex = GetNParameters();

   // add a parameter for the expectation value
   AddParameter(name, Nmin, Nmax);

   // get efficiency parameter index
   int effindex = GetNParameters();

    // add a parameter for the efficiency
   AddParameter(Form("Efficiency_%s", name), -5.0, 5.0);

   // add histogram, name and index to containers
   fTemplateHistogramContainer.push_back(hist);
   fHistogramIndexContainer.push_back(GetNTemplatesType(0));
   fFunctionIndexContainer.push_back(-1);
   fTemplateTypeContainer.push_back(0);
   fTemplateNameContainer.push_back(name);
   fTemplateParIndexContainer.push_back(parindex);
   fEffParIndexContainer.push_back(effindex);

   // set efficiency histograms to one without uncertainty
   fEffHistogramContainer.push_back(TH1D());
   fEffErrHistogramContainer.push_back(TH1D());
   SetTemplateEfficiency(name, 1., 0.);

   // set prior for efficiencies
   SetPriorGauss(Form("Efficiency_%s", name), 0., 1.);

   // add systematic uncertainties
   for (int i = 0; i < GetNSystErrors(); ++i) {
      std::vector<TH1D> histvector = fSystErrorHistogramContainer.at(i);
      histvector.push_back(histsysterror);
   }

   // successfully added histogram to container
   return 1;
}

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

Definition at line 350 of file BCTemplateFitter.cxx.

{
   // check if prior makes sense
   if (fFlagPhysicalLimits && Nmin < 0)
      Nmin = 0;

   if (Nmin > Nmax) {
      BCLog::OutError("BCTemplateFitter::AddTemplate : Lower limit exceeds upper limit.");
      return 0;
   }

   // get parameter index
   int parindex = GetNParameters();

   // add a parameter for the expectation value
   AddParameter(name, Nmin, Nmax);

   // get efficiency parameter index
   int effindex = GetNParameters();

    // add a parameter for the efficiency
   AddParameter(Form("Efficiency_%s", name), -5.0, 5.0);

   // add all additional parameters
   int npar = func.GetNpar();

   for (int i = 0; i < npar; ++i) {
      double parmin;
      double parmax;
      func.GetParLimits(i, parmin, parmax);
      AddParameter(func.GetParName(i), parmin, parmax);
   }


   // add histogram, name and index to containers
   fTemplateFunctionContainer.push_back(func);
   fHistogramIndexContainer.push_back(-1);
   fFunctionIndexContainer.push_back(GetNTemplatesType(1));
   fTemplateTypeContainer.push_back(1);
   fTemplateNameContainer.push_back(name);
   fTemplateParIndexContainer.push_back(parindex);
   fEffParIndexContainer.push_back(effindex);

   // set efficiency histograms to one without uncertainty
   fEffHistogramContainer.push_back(TH1D());
   fEffErrHistogramContainer.push_back(TH1D());
   SetTemplateEfficiency(name, 1., 0.);

   // set prior for efficiencies
   SetPriorGauss(Form("Efficiency_%s", name), 0., 1.);

   // successfully added histogram to container
   return 1;
}

double BCTemplateFitter::CalculateChi2

(

const std::vector< double > &

parameters

)

Calculates and returns the chi2 value. The chi2 is calculated using the expectation value for the uncertainty.

Definition at line 975 of file BCTemplateFitter.cxx.

{
   double chi2 = 0;

   // loop over bins
   for (int ibin = 1; ibin <= fNBins; ++ibin) {
      // get expectation
      double nexp = Expectation(ibin, parameters);

      // get data
      double ndata = fHistData.GetBinContent(ibin);

      // add to chi2
      chi2 += (nexp - ndata) * (nexp - ndata) / nexp;
   }

   // return chi2
   return chi2;
}

double BCTemplateFitter::CalculateChi2Prob

(

const std::vector< double > &

parameters

)

Calculates and returns the chi2-probability.

Definition at line 996 of file BCTemplateFitter.cxx.

{
   double chi2 = CalculateChi2(parameters);
   int ndf = GetNDF();

   // return chi2 probability
   return TMath::Prob(chi2, ndf);
}

double BCTemplateFitter::CalculateKSProb

(

)

Calculates and returns the Kolmogorov-Smirnov-probability.

Definition at line 1013 of file BCTemplateFitter.cxx.

{
   // create a temporary histogram with the sum of all contributions
    TH1D * histsum = new TH1D(fHistData);

   // get best fit parameters
   std::vector<double> parameters = GetBestFitParameters();

   // loop over bins
   for (int ibin = 1; ibin <= fNBins; ++ibin) {
      // get expectation
      double nexp = Expectation(ibin, parameters);

      // set bin content
      histsum->SetBinContent(ibin, nexp);
   }

   // perform KS test
   double ksprob = histsum->KolmogorovTest(&fHistData);

   // delete histogram
   delete histsum;

   return ksprob;
}

double BCTemplateFitter::CalculateMaxLike

(

)

Calculates and returns the Likelihood at the global mode.

Definition at line 1006 of file BCTemplateFitter.cxx.

{
   // return maximum likelihood
   return Eval( GetBestFitParameters() );
}

double BCTemplateFitter::CalculatePValue

(

)

Calculates and returns the p-value for the global mode.

Definition at line 1040 of file BCTemplateFitter.cxx.

{
   // get best fit parameters
   std::vector<double> par = GetBestFitParameters();

   // check size of parameter vector
   if (par.size() != GetNParameters()) {
      BCLog::OutWarning("BCBCTemplateFitter::CalculatePValueFast : Number of parameters is inconsistent.");
      return -1;
   }

   std::vector<int> histogram;
   std::vector<double> expectation;
   histogram.assign(fNBins, 0);
   expectation.assign(fNBins, 0);

   double logp = 0;
   double logp_start = 0;
   int counter_pvalue = 0;

   // define starting distribution
   for (int ibin = 1; ibin <= fNBins; ++ibin) {

      // get expectation
      double nexp = Expectation(ibin, par);

      /*
      double nexp = 0;

      // loop over all templates
      for (int itemp = 0; itemp < ntemplates; ++itemp) {
         int tempindex = fTemplateParIndexContainer.at(itemp);
         int effindex = fEffParIndexContainer.at(itemp);

         // get efficiency for the bin
         double efficiency = fEffHistogramContainer.at(itemp).GetBinContent(ibin);

         // modify efficiency by uncertainty
         double efferr = fEffErrHistogramContainer.at(itemp).GetBinContent(ibin);

         // check efficiency error
         if (efferr > 0)
            efficiency = TMath::Max(0., efficiency + par.at(effindex) * efferr);

         // add expectation from bin
         nexp += par.at(tempindex) * efficiency * fTemplateHistogramContainer.at(itemp).GetBinContent(ibin);
      }
      */


      histogram[ibin-1]   = int(nexp);
      expectation[ibin-1] = nexp;

      // calculate p;
      logp += BCMath::LogPoisson(double(int(nexp)), nexp);
      logp_start += BCMath::LogPoisson(fHistData.GetBinContent(ibin), nexp);
   }

   int niter = 100000;

   // loop over iterations
   for (int iiter = 0; iiter < niter; ++iiter) {
      // loop over bins
      for (int ibin = 0; ibin < fNBins; ++ibin) {
         // random step up or down in statistics for this bin
         double ptest = fRandom->Rndm() - 0.5;

         // increase statistics by 1
         if (ptest > 0) {
            // calculate factor of probability
            double r = expectation.at(ibin) / double(histogram.at(ibin) + 1);

            // walk, or don't (this is the Metropolis part)
            if (fRandom->Rndm() < r) {
               histogram[ibin] = histogram.at(ibin) + 1;
               logp += log(r);
            }
         }

         // decrease statistics by 1
         else if (ptest <= 0 && histogram[ibin] > 0) {
            // calculate factor of probability
            double r = double(histogram.at(ibin)) / expectation.at(ibin);

            // walk, or don't (this is the Metropolis part)
            if (fRandom->Rndm() < r) {
               histogram[ibin] = histogram.at(ibin) - 1;
               logp += log(r);
            }
         }
      } // end of looping over bins

      // increase counter
      if (logp < logp_start)
         counter_pvalue++;

   } // end of looping over iterations

   // calculate p-value
   return double(counter_pvalue) / double(niter);
}

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

Add the calculation of a certain ratio.

Parameters:

	index	Index in the numerator.
	indices	Vector of indices in the denominator.
	rmin	The minimum ratio
	rmax	The maximum ratio

Returns:

An error code.

Definition at line 550 of file BCTemplateFitter.cxx.

{
   // get number of templates
   int ntemplates = GetNTemplates();

   // check index
   if (index < 0 || index >= ntemplates) {
      return 0;
   }

   // check indices
   for (int i = 0; i < int(indices.size()); ++i) {
      if (indices.at(i) < 0 || indices.at(i) >= ntemplates) {
         return 0;
      }
   }

   // create temporary vector
   std::vector<int> tempvec;
   tempvec.push_back(index);
   for (int i = 0; i < int(indices.size()); ++i)
      tempvec.push_back(indices.at(i));

   // add ratio
   fIndicesRatios1D.push_back(tempvec);

   // get number of ratios
   int nratios = int(fHistRatios1D.size());

   // create histogram
   double fmin = rmin;
   double fmax = rmax;
   if (fFlagPhysicalLimits) {
      fmin = TMath::Max(rmin, 0.0);
      fmax = TMath::Min(1.0, rmax);
   }

   TH1D hist_ratio1d(TString::Format("ratio %i", nratios), ";r;p(r|data)", 100, fmin, fmax);
   fHistRatios1D.push_back(hist_ratio1d);

   // no error
   return 1;
}

TH1D BCTemplateFitter::CombineUncertainties

(

const char *

name

)

Cobine all sources of systematic uncertainties (including efficiency) by summing the squares.

Parameters:

name

The name of the template

Returns:

A histogram with the total uncertainty

Definition at line 1454 of file BCTemplateFitter.cxx.

{
   // get number of sources of systematic uncertainty
   int nsyst = GetNSystErrors();

   // get template index
   int tempindex = GetIndexTemplate(name);

   // create new histogram
   //    TH1D hist = TH1D("", "", fNBins, fXmin, fXmax);
    TH1D hist = TH1D(fHistData);

   // fill histogram
   for (int ibin = 1; ibin <= fNBins; ++ibin) {

      // define total uncertainty
      double err = 0;

      // add efficiency uncertainty squared
      double erreff = fEffErrHistogramContainer.at(tempindex).GetBinContent(ibin);
      err += erreff * erreff;

      // add systematic uncertainty squared
      for (int isyst = 0; isyst < nsyst; ++isyst) {
         double errsyst = fSystErrorHistogramContainer.at(isyst).at(tempindex).GetBinContent(ibin);
         err += errsyst*errsyst;
      }

      // take square root
      err = sqrt(err);

      // set bin content
      hist.SetBinContent(ibin, err);
   }

   // return histogram
   return hist;
}

int BCTemplateFitter::CompareHistogramProperties	(	TH1D	hist1,
		TH1D	hist2
	)

Checks if two histograms have the same properties.

Returns:

0 if not, 1 if they have the same properties.

Definition at line 595 of file BCTemplateFitter.cxx.

{
   // compare number of bins
   if (hist1.GetNbinsX() != hist2.GetNbinsX())
      return 0;

   // compare minimum x-values
   if (hist1.GetXaxis()->GetXmin() != hist2.GetXaxis()->GetXmin())
      return 0;

   // compare maximum x-values
   if (hist1.GetXaxis()->GetXmax() != hist2.GetXaxis()->GetXmax())
      return 0;

   // conclusion: they have the same properties
   return 1;
}

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

Add a correlation among two sources of systematic uncertainty.

Parameters:

	errorname1	The name of the first source.
	errorname2	The name of the second source.
	corr	The correlation coefficiency.

Returns:

An error code. Constrains a sum of contributions. Assume a Gaussian prior.

Parameters:

	indices	The vector of indicies for the contributions.
	mean	The mean value of the prior.
	rms	The standard deviation of the prior.

Returns:

An error code.

Definition at line 1276 of file BCTemplateFitter.cxx.

{
   // add contraint to container(s)
   fConstraintSumIndices.push_back(indices);
   fConstraintSumMean.push_back(mean);
   fConstraintSumRMS.push_back(rms);

   // no error
   return 1;
}

TH1D BCTemplateFitter::CreateErrorHist	(	TH1D	hist,
		TH1D	histerr
	)

Create a histogram with specified uncertainties.

Parameters:

	hist	The histogram to be copied.
	histerr	The uncertainties of the new histogram.

Returns:

A histogram with uncertainties.

Definition at line 1426 of file BCTemplateFitter.cxx.

{
   // check histogram properties
   if (CompareHistogramProperties(fHistData, hist) != 1) {
      BCLog::OutError("BCTemplateFitter::CreateErrorHist : Histograms are incompatible.");
      return hist;
   }

   // copy histogram
   TH1D h = hist;

   // set style
   h.SetStats(kFALSE);
   h.SetFillColor(kYellow);
   h.SetFillStyle(1001);

   // get number of bins
   int nbins = hist.GetNbinsX();

   // loop over bins
   for (int i = 1; i <= nbins; ++i)
      h.SetBinError(i, histerr.GetBinContent(i) * hist.GetBinContent(i));

   // return histogram
   return h;
}

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

Return the expectation value in a certain bin.

Parameters:

	binnumber	The number of the bin.
	parameters	The parameter values.

Definition at line 98 of file BCTemplateFitter.cxx.

{
   // initialize number of expected events with 0
   double nexp = 0;

   // get number of templates
   int ntemplates = GetNTemplates();

   // loop over all templates
   for (int itemp = 0; itemp < ntemplates; ++itemp) {

      // get template index
      int templateindex = fTemplateParIndexContainer.at(itemp);

      // calculate efficiency for this template and bin
      double efficiency = TemplateEfficiency(itemp, binnumber, parameters);

      // calculate probability for this template and bin
      double probability = TemplateProbability(itemp, binnumber, parameters);

      // calculate expectation value for this template
      nexp += parameters.at(templateindex)
         * efficiency
         * probability;
   }

   // check that expectation is larger or equal to zero
   if (nexp < 0) {
      BCLog::OutWarning("BCTemplateFitter::Expectation : Expectation value smaller than 0. Force it to be 0.");
      nexp = 0;
   }

   return nexp;
}

int BCTemplateFitter::FixTemplateFunctions

(

const std::vector< double > &

parameters

)

Definition at line 200 of file BCTemplateFitter.cxx.

{
   // get number of templates
   int ntemplates = GetNTemplates();

   // fix all function parameters
   for (int i = 0; i < ntemplates; ++i) {
      if (fTemplateTypeContainer.at(i) == 1) {
         int index = GetFunctionIndex(i);
         TF1* func = &(fTemplateFunctionContainer.at(index));
         int npar = func->GetNpar();
         int tempindex = GetParIndexTemplate(i);
         for (int i = 0; i < npar; ++i) {
            func->FixParameter(i, parameters.at(tempindex+2+i));
         }
      }
   }

   // no errors
   return 1;
}

TH1D BCTemplateFitter::GetData

(

)

[inline]

Return a template histogram.

Parameters:

name

The template name. Return the histogram containing the data.

Definition at line 164 of file BCTemplateFitter.h.

         { return fHistData; }

int BCTemplateFitter::GetFunctionIndex

(

int

templatenumber

)

[inline]

Get function index of template. The number of the template

Definition at line 178 of file BCTemplateFitter.h.

         { return fFunctionIndexContainer.at(templatenumber); };

int BCTemplateFitter::GetHistogramIndex

(

int

templatenumber

)

[inline]

Get histogram index of template. The number of the template

Definition at line 171 of file BCTemplateFitter.h.

         { return fHistogramIndexContainer.at(templatenumber); };

TH1D BCTemplateFitter::GetHistRatio1D

(

int

index

)

[inline]

Return a ratio histogram

Definition at line 110 of file BCTemplateFitter.h.

         { return fHistRatios1D.at(index); }

std::vector<TH1D> BCTemplateFitter::GetHistRatios1D

(

)

[inline]

Return the vector of histgrams containing the ratios.

Definition at line 92 of file BCTemplateFitter.h.

         { return fHistRatios1D; }

int BCTemplateFitter::GetIndexSystError

(

const char *

name

)

Return the index of a source of systematic uncertainty.

Parameters:

name

The name of the source.

Definition at line 1513 of file BCTemplateFitter.cxx.

{
   int index = -1;
   int n = GetNSystErrors();

   for (int i = 0; i < n; i++)
      if (name == fSystErrorNameContainer.at(i))
         index = i;

   if (index < 0) {
      BCLog::OutWarning("BCTemplateFitter::GetIndexSystError : Template does not exist.");
      return 0;
   }

   // return index
   return index;
}

int BCTemplateFitter::GetIndexTemplate

(

const char *

name

)

Return the index of a template.

Parameters:

name

The template name.

Definition at line 1494 of file BCTemplateFitter.cxx.

{
   int index = -1;
   int n = GetNTemplates();

   for (int i = 0; i < n; i++)
      if (name == fTemplateNameContainer.at(i))
         index = i;

   if (index < 0) {
      BCLog::OutWarning("BCTemplateFitter::GetIndexTemplate : Template does not exist.");
      return 0;
   }

   // return index
   return index;
}

int BCTemplateFitter::GetNDF

(

)

[inline]

Return the number of degrees of freedom.

Definition at line 80 of file BCTemplateFitter.h.

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

int BCTemplateFitter::GetNRatios

(

)

[inline]

Return the number of ratios to calculate.

Definition at line 86 of file BCTemplateFitter.h.

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

int BCTemplateFitter::GetNSystErrors

(

)

[inline]

Return the number of sources of systematic uncertainties.

Definition at line 74 of file BCTemplateFitter.h.

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

int BCTemplateFitter::GetNTemplates

(

)

[inline]

Return the number of templates.

Definition at line 62 of file BCTemplateFitter.h.

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

int BCTemplateFitter::GetNTemplatesType

(

int

templatetype

)

Return the number of templates of a certain type

Parameters:

templatetype

the type of the template

Definition at line 1623 of file BCTemplateFitter.cxx.

{
   // get number of templates
   int ntemplates = GetNTemplates();

   // initialize counter to zero
   int counter = 0;

   // loop over all templates
   for (int i = 0; i < ntemplates; ++i) {
      // increase counter if template of same type
      if (templatetype == fTemplateTypeContainer.at(i))
         counter++;
   }

   // return counter
   return counter;
}

int BCTemplateFitter::GetParIndexEff

(

const char *

name

)

Return the parameter index corresponding to an efficiency.

Parameters:

name

The name of the template associated with the efficiency.

Definition at line 1566 of file BCTemplateFitter.cxx.

{
   int index = -1;
   int n = GetNTemplates();

   for (int i = 0; i < n; i++)
      if (name == fTemplateNameContainer.at(i))
         index = fTemplateParIndexContainer.at(i);

   if (index < 0) {
      BCLog::OutWarning("BCTemplateFitter::GetParIndexEff : Template does not exist.");
      return 0;
   }

   // return index
   return index;
}

int BCTemplateFitter::GetParIndexSystError

(

const char *

name

)

Return the parameter index corresponding to an efficiency.

Parameters:

index

The index of the template associated with the efficiency.

Definition at line 1585 of file BCTemplateFitter.cxx.

{
   int index = -1;
   int n = GetNTemplates();

   for (int i = 0; i < n; i++)
      if (name == fSystErrorNameContainer.at(i))
         index = fSystErrorParIndexContainer.at(i);

   if (index < 0) {
      BCLog::OutWarning("BCTemplateFitter::GetParIndexStatError : Systematic error does not exist.");
      return 0;
   }

   // return index
   return index;
}

int BCTemplateFitter::GetParIndexTemplate

(

const char *

name

)

Return the parameter index corresponding to a template.

Parameters:

name

The template name.

Definition at line 1532 of file BCTemplateFitter.cxx.

{
   int index = -1;
   int n = GetNTemplates();

   for (int i = 0; i < n; i++)
      if (name == fTemplateNameContainer.at(i))
         index = fTemplateParIndexContainer.at(i);

   if (index < 0) {
      BCLog::OutWarning("BCTemplateFitter::GetParIndexTemplate : Template does not exist.");
      return 0;
   }

   // return index
   return index;
}

int BCTemplateFitter::GetParIndexTemplate

(

int

index

)

Return the parameter index corresponding to a template.

Parameters:

index

The template index.

Definition at line 1551 of file BCTemplateFitter.cxx.

{
   // get number of templates
   int n = GetNTemplates();

   if (index < 0 || index > n) {
      BCLog::OutError("BCTemplateFitter::GetParIndexTemplate : Index out of range.");
      return -1;
   }

   // return index
   return fTemplateParIndexContainer.at(index);
}

std::vector<TF1> BCTemplateFitter::GetTemplateFunctionContainer

(

)

[inline]

Return the container of template functions.

Definition at line 104 of file BCTemplateFitter.h.

         { return fTemplateFunctionContainer; };

std::vector<TH1D> BCTemplateFitter::GetTemplateHistogramContainer

(

)

[inline]

Return the container of template histograms.

Definition at line 98 of file BCTemplateFitter.h.

         { return fTemplateHistogramContainer; };

std::vector<int> BCTemplateFitter::GetTemplateParIndexContainer

(

)

[inline]

Return container of parameter indeces for templates

Definition at line 184 of file BCTemplateFitter.h.

         { return fTemplateParIndexContainer; };

TH2D* BCTemplateFitter::GetUncertaintyHistogramExp

(

)

[inline]

Definition at line 113 of file BCTemplateFitter.h.

      { return fUncertaintyHistogramExp; };

TH2D* BCTemplateFitter::GetUncertaintyHistogramObsPosterior

(

)

[inline]

Definition at line 116 of file BCTemplateFitter.h.

      { return fUncertaintyHistogramObsPosterior; };

int BCTemplateFitter::Initialize

(

)

Initialize the fitting procedure.

Returns:

An error code.

Definition at line 490 of file BCTemplateFitter.cxx.

{
   // check data integral
   if (fHistData.Integral() <= 0) {
      BCLog::OutError("BCTemplateFitter::Initialize : Normalization of data histogram is zero or less than that.");
      return 0;
   }

   // create histograms for uncertainty determination
   //   double maximum = 1.5 * fHistData.GetMaximum();
   double minimum = TMath::Max(0., fHistData.GetMinimum() - 5.*sqrt(fHistData.GetMinimum()));
   double maximum = fHistData.GetMaximum() + 5.*sqrt(fHistData.GetMaximum());

   std::vector<double> a(fHistData.GetNbinsX()+1);
   for (int i = 0; i < fHistData.GetNbinsX()+1; ++i) {
      a[i] = fHistData.GetXaxis()->GetBinLowEdge(i+1);
   }

   fUncertaintyHistogramExp = new TH2D(
         TString::Format("UncertaintyExp_%i", BCLog::GetHIndex()), "",
         //         fHistData.GetNbinsX(), fHistData.GetXaxis()->GetXbins()->GetArray(),
         fHistData.GetNbinsX(), &a[0],
         //         fHistData.GetNbinsX(), fHistData.GetXaxis()->GetXmin(), fHistData.GetXaxis()->GetXmax(),
         100, minimum, maximum);

   fUncertaintyHistogramObsPosterior = new TH2D(
         TString::Format("UncertaintyObsPosterior_%i", BCLog::GetHIndex()), "",
         fHistData.GetNbinsX(), &a[0],
         //         fHistData.GetNbinsX(), fHistData.GetXaxis()->GetXbins()->GetArray(),         //         fHistData.GetNbinsX(), fHistData.GetXaxis()->GetXmin(), fHistData.GetXaxis()->GetXmax(),
         int(maximum) + 1, -0.5, double(int(maximum))+0.5);

   // create histogram containing the normalization
   double xmin = 0;
   double xmax = 0;
   int ntemplates = GetNTemplates();

   // calculate the limits on the norm from the sum of all parameter
   // limits
   for (int i = 0; i < ntemplates; ++i) {
      // get parameter index
      int parindex = fTemplateParIndexContainer.at(i);
      int effindex = fEffParIndexContainer.at(i);

      // get parameter
      BCParameter * par = this->GetParameter(parindex);
      BCParameter * eff = this->GetParameter(effindex);

      // increate limits
      xmin += par->GetLowerLimit() * eff->GetLowerLimit();
      xmax += par->GetUpperLimit() * eff->GetUpperLimit();
   }

   // create new histogram for norm
   fHistNorm = TH1D("", ";N_{norm};dN/dN_{norm}", 100, xmin, xmax);

   // no error
   return 1;
}

double BCTemplateFitter::LogLikelihood

(

const std::vector< double > &

parameters

)

[virtual]

Calculates and returns the log of the Likelihood at a given point in parameter space.

Implements BCModel.

Definition at line 74 of file BCTemplateFitter.cxx.

{
   double logprob = 0.;

   // fix the parameters of all functions
   FixTemplateFunctions(parameters);

   // loop over bins
   for (int ibin = 1; ibin <= fNBins; ++ibin) {
      // get expectation
      double nexp = Expectation(ibin, parameters);

      // get data
      double ndata = fHistData.GetBinContent(ibin);

      // add Poisson term
      logprob += BCMath::LogPoisson(ndata, nexp);
   }

   // return log likelihood
   return logprob;
}

void BCTemplateFitter::MCMCUserIterationInterface

(

)

[virtual]

Overloaded from BCIntegrate.

Reimplemented from BCIntegrate.

Definition at line 1143 of file BCTemplateFitter.cxx.

{
   // loop over all bins
   for (int ibin = 1; ibin <= fNBins; ++ibin) {
      // get bin content
      double bincontent = Expectation(ibin, fMCMCx);

      // set bin content
      fUncertaintyHistogramExp->Fill(fHistData.GetBinCenter(ibin), bincontent);

      // loop over bins in the other direction
      int nbinsy = fUncertaintyHistogramObsPosterior->GetNbinsY();
      for (int jbin = 1; jbin <= nbinsy; ++jbin) {
         int n = jbin - 1;
         if (fabs(n - bincontent) < 2*sqrt(bincontent))
            fUncertaintyHistogramObsPosterior->Fill(fHistData.GetBinCenter(ibin), n, TMath::Poisson(bincontent, n));
      }
   }

   // fill normalization
   fHistNorm.Fill(fNorm);

   // fill ratios
   int nratios = int( fIndicesRatios1D.size() );

   // loop over fractions to fill
   for (int i = 0; i < nratios; ++i) {
      int nsum = int( (fIndicesRatios1D.at(i)).size() ) - 1;
      double sum = 0;
      for (int j = 1; j <= nsum; ++j) {
         int indexsum = fIndicesRatios1D.at(i).at(j);
         sum += fMCMCx.at(fTemplateParIndexContainer.at(indexsum));
      }

      fHistRatios1D.at(i).Fill(fMCMCx.at(fTemplateParIndexContainer.at(fIndicesRatios1D.at(i).at(0)))/sum);
   }

}

int BCTemplateFitter::PerformFit

(

)

Perform the template fit.

Returns:

An error code.

Definition at line 1604 of file BCTemplateFitter.cxx.

{
   // initialize
   if (!Initialize()) {
      BCLog::OutError("BCTemplateFitter::PerformFit : Could not initialize template fitter.");
      return 0;
   }

   // run Markov Chains
   MarginalizeAll();

   // find global mode
   FindMode();

   // no error
   return 1;
}

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

Print the ratios and the norm.

Parameters:

	filename	The filename.
	option	Plot options
	pvalue	Value which goes with plot options (see BAT manual).

Definition at line 1183 of file BCTemplateFitter.cxx.

{
   int nratios = int(fHistRatios1D.size());

   TCanvas * c1 = new TCanvas("");

   TPostScript * ps = new TPostScript(filename, 112);
   ps->NewPage();

   c1->cd();
   BCH1D * h1temp = new BCH1D(&fHistNorm);
   h1temp->Draw();
   c1->Update();
   ps->NewPage();
   for (int i = 0; i < nratios; ++i) {
      c1->Update();
      ps->NewPage();
      c1->cd();
      BCH1D* h1temp = new BCH1D(&fHistRatios1D.at(i));
      h1temp->Draw(options, ovalue);
   }
   c1->Update();
   ps->Close();

   // free memory
   delete c1;
   delete ps;
}

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

Prints the stack of templates scaled with the global mode. The data is plotted on top. The following options are available:

"L" : adds a legend

"E0" : symmetric errorbars on the data points with a length of sqrt(obs).

"E1" : symmetric errorbars on the sum of templates with a length of sqrt(exp).

"E2" : asymmetric errorbars on the sum of templates from error propagation of the parameters.

"E3" : asymmetric errorbars on the data points. The errorbars mark the 16% and 85% quantiles of the Poisson distribution rounded to the lower or upper integer, respectively.

Parameters:

	filename	The name of the file the output is printed to.
	options	Plotting options

Definition at line 614 of file BCTemplateFitter.cxx.

{
   int nbins = fHistData.GetNbinsX();
   int ntemplates = GetNTemplates();

   // check options
   bool flag_legend = false;
   bool flag_error0 = false; // symm. poisson error for data
   bool flag_error1 = false; // symm. poisson error for exp.
   bool flag_error2 = false; // asymm. poisson error of expectation value
   bool flag_error3 = false; // asymm. poisson error of expected no. of events
   bool flag_diff   = false; // plot difference between data and expectation below stack plot
   bool flag_logx   = false; // plot x-axis in log-scale
   bool flag_logy   = false; // plot y-axis in log-scale

   if (std::string(options).find("L") < std::string(options).size())
      flag_legend = true;

   if (std::string(options).find("E0") < std::string(options).size())
      flag_error0 = true;

   if (std::string(options).find("E1") < std::string(options).size())
      flag_error1 = true;

   if (std::string(options).find("E2") < std::string(options).size())
      flag_error2 = true;

   if (std::string(options).find("E3") < std::string(options).size())
      flag_error3 = true;

   if (std::string(options).find("D") < std::string(options).size())
      flag_diff = true;

   if (std::string(options).find("logx") < std::string(options).size())
      flag_logx = true;

   if (std::string(options).find("logy") < std::string(options).size())
      flag_logy = true;

   // create canvas
   TCanvas* c1 = new TCanvas("", "", 700, 700);
   c1->cd();
   TPad * pad1;
   TPad * pad2;

   double fraction_pads = 0.3;

   if(!flag_diff)
      fraction_pads=0.0;

   if (flag_diff) {
      pad1 = new TPad("pad1", "", 0.0, fraction_pads, 1.0, 1.0);
      pad1->SetTopMargin   (0.13/(1.0-fraction_pads));
      pad1->SetBottomMargin(0.0);
      pad1->SetLeftMargin  (0.15);
      pad1->SetRightMargin (0.13);
      pad1->SetFillColor(kWhite);
      pad2 = new TPad("pad2", "", 0.0, 0.0, 1.0, fraction_pads);
      pad2->SetTopMargin   (0.0);
      pad2->SetBottomMargin(0.15 / fraction_pads);
      pad2->SetLeftMargin  (0.15);
      pad2->SetRightMargin (0.13);
      pad2->SetFillColor(kWhite);
      if (flag_logx) {
         pad1->SetLogx();
         pad2->SetLogx();
      }
      if (flag_logy)
         pad1->SetLogy();
      pad1->Draw();
      pad2->Draw();
   }
   else {
      pad1 = new TPad("pad1", "",0.0, 0.0, 1.0, 1.0);
      pad1->SetFillColor(kWhite);
      pad2 = new TPad();
      if (flag_logx)
         pad1->SetLogx();
      if (flag_logy)
         pad1->SetLogy();
      pad1->Draw();
   }

   pad1->cd();

   // set style and draw data
   double ymin = 0.01;
   double ymax = 1.1 * (fHistData.GetMaximum() + sqrt(fHistData.GetMaximum()));
   fHistData.GetYaxis()->SetRangeUser(ymin, ymax);
   fHistData.GetXaxis()->SetNdivisions(505);
   if (flag_diff) {
      fHistData.GetXaxis()->SetLabelSize(fHistData.GetXaxis()->GetLabelSize()/(1.0-fraction_pads));
      fHistData.GetXaxis()->SetLabelOffset(fHistData.GetXaxis()->GetLabelOffset()*(1.0-fraction_pads));
      fHistData.GetXaxis()->SetTitleSize(fHistData.GetXaxis()->GetTitleSize()/(1.0-fraction_pads));
      fHistData.GetXaxis()->SetTitleOffset(fHistData.GetXaxis()->GetTitleOffset()*(1.0-fraction_pads));
      fHistData.GetYaxis()->SetLabelSize(fHistData.GetYaxis()->GetLabelSize()/(1.0-fraction_pads));
      fHistData.GetYaxis()->SetLabelOffset(fHistData.GetYaxis()->GetLabelOffset()/(fraction_pads));
      fHistData.GetYaxis()->SetTitleSize(fHistData.GetYaxis()->GetTitleSize()/(1.0-fraction_pads));
      fHistData.GetYaxis()->SetTitleOffset(fHistData.GetYaxis()->GetTitleOffset()*(1.0-fraction_pads));
   }
   fHistData.Draw("P");

   // create a histogram with the sum of all contributions
   //    TH1D * histsum = (TH1D*) fHistData.Clone("temp");
    TH1D * histsum = new TH1D(fHistData);

   // create stack
   THStack stack("histostack","");

   // create legends
   TLegend* legend1;
   TLegend* legend2;

   if (flag_diff)
      legend1 = new TLegend(0.15, (0.88-fraction_pads)/(1-fraction_pads), 0.50, 0.99);
   else
      legend1 = new TLegend(0.15, 0.88, 0.50, 0.99);
   legend1->SetBorderSize(0);
   legend1->SetFillColor(kWhite);
   legend1->AddEntry(&fHistData, "Data", "LEP");
   legend1->AddEntry(&fHistData, "Total expected uncertainty", "LE");

   double y = 0.99;
   if (ntemplates > 2 && ntemplates <7)
      y -= 0.11 / 4. * double(ntemplates - 2);
   legend2 = new TLegend(0.50,(y-fraction_pads)/(1-fraction_pads) , 0.85, 0.99);
   legend2->SetBorderSize(0);
   legend2->SetFillColor(kWhite);


   // create new histograms for plotting
   std::vector<TH1D*> histcontainer;

   // get best fit parameters
   std::vector<double> parameters = GetBestFitParameters();

   // fix the parameters of all functions
   FixTemplateFunctions(parameters);

   // scale histograms and add to stack and legend
   for (int itemp = 0; itemp < ntemplates; ++itemp) {

      // get template index
      int templateindex = fTemplateParIndexContainer.at(itemp);

      // create new histogram
      TH1D* hist;
      if (fTemplateTypeContainer[itemp] == 0) {
         int histogramindex = GetHistogramIndex(itemp);
         hist = new TH1D( fTemplateHistogramContainer.at(histogramindex) );
      }
      else if (fTemplateTypeContainer[itemp] == 1) {
         hist = new TH1D( fHistData );
         for (int i = 1; i <= fNBins; ++i) {
            double bincenter = fHistData.GetBinCenter(i);
            int index = GetFunctionIndex(itemp);
            double bincontent = fTemplateFunctionContainer.at(index).Eval(bincenter);
            hist->SetBinContent(i, bincontent);
         }
         // scale histogram
         hist->Scale(1.0/hist->Integral());
      }

      // set histogram color and style
      hist->SetFillColor(2 + itemp);
      hist->SetFillStyle(1001);
      hist->SetStats(kFALSE);

      // add histogram to list of histograms
      histcontainer.push_back(hist);

      // loop over bins and set them
      for (int ibin = 1; ibin <= fNBins; ++ibin) {

         // calculate efficiency for this template and bin
         double efficiency = TemplateEfficiency(itemp, ibin, parameters);

         // calculate probability for this template and bin
         double probability = TemplateProbability(itemp, ibin, parameters);

         // calculate expectation value for this template
         double nexp = parameters.at(templateindex)
            * efficiency
            * probability;

         // set bin content
         hist->SetBinContent(ibin, nexp);
      }


      // add histogram to stack
      stack.Add(hist);
      if (itemp < 2)
         legend1->AddEntry(hist, fTemplateNameContainer.at(itemp).data(), "F");
      else if (itemp < 6)
         legend2->AddEntry(hist, fTemplateNameContainer.at(itemp).data(), "F");
   }

   // loop over all bins
   for (int ibin = 1; ibin <= nbins; ++ibin) {
      // set bin content
      histsum->SetBinContent(ibin, Expectation(ibin, parameters));
   }

   // define error graph
   TGraphAsymmErrors * graph_error_exp = new TGraphAsymmErrors(nbins);
   graph_error_exp->SetLineWidth(2);
   graph_error_exp->SetMarkerStyle(0);

   TGraphAsymmErrors * graph_error_obs = new TGraphAsymmErrors(nbins);
   graph_error_obs->SetMarkerStyle(0);

   // calculate uncertainty
   if (flag_error1)
      for (int i = 1; i <= nbins; ++i) {
         double nexp = histsum->GetBinContent(i);
         histsum->SetBinError(i, sqrt(nexp));
         histsum->SetMarkerStyle(0);
      }

   if (flag_error2)
      for (int i = 1; i <= nbins; ++i) {
         TH1D * proj = fUncertaintyHistogramExp->ProjectionY("_py", i, i);
         if (proj->Integral() > 0)
            proj->Scale(1.0 / proj->Integral());
         double quantiles[3];
         double sums[3] = {0.16, 0.5, 0.84};
         proj->GetQuantiles(3, quantiles, sums);
         graph_error_exp->SetPoint(i-1, histsum->GetBinCenter(i), quantiles[1]);
         graph_error_exp->SetPointError(i-1, 0.0, 0.0, quantiles[1] - quantiles[0], quantiles[2]-quantiles[1]);
         delete proj;
      }

   if (flag_error3)
      for (int i = 1; i <= nbins; ++i) {
         TH1D * proj = fUncertaintyHistogramObsPosterior->ProjectionY("_py", i, i);
         if (proj->Integral() > 0)
            proj->Scale(1.0 / proj->Integral());
         double quantiles[3];
         double sums[3] = {0.16, 0.5, 0.84};
         proj->GetQuantiles(3, quantiles, sums);
         graph_error_obs->SetPoint(i-1, histsum->GetBinCenter(i), quantiles[1]);
         graph_error_obs->SetPointError(i-1, 0.0, 0.0, quantiles[1] - TMath::Floor(quantiles[0]), TMath::Ceil(quantiles[2])-quantiles[1]);
         delete proj;
      }

   // create difference histogram
   TH1D * hist_diff = 0;

   TGraphAsymmErrors * graph_diff_exp = 0;

   if (flag_diff) {
      ymin = 0;
      ymax = 0;
      //      hist_diff = new TH1D("hist_diff", "", nbins, histsum->GetXaxis()->GetXmin(), histsum->GetXaxis()->GetXmax() );
      //      hist_diff = new TH1D(*histsum);
      std::vector<double> a(fHistData.GetNbinsX()+1);
      for (int i = 0; i < fHistData.GetNbinsX()+1; ++i) {
         a[i] = fHistData.GetXaxis()->GetBinLowEdge(i+1);
      }

      hist_diff = new TH1D("hist_diff", "", nbins, &a[0]);
      hist_diff->SetName("hist_diff");
      hist_diff->GetXaxis()->SetTitle(fHistData.GetXaxis()->GetTitle());
      hist_diff->GetYaxis()->SetTitle("#Delta N");
      hist_diff->GetXaxis()->SetNdivisions(505);
      hist_diff->GetXaxis()->SetLabelSize(hist_diff->GetXaxis()->GetLabelSize()/(fraction_pads));
      hist_diff->GetXaxis()->SetLabelOffset(hist_diff->GetXaxis()->GetLabelOffset()/fraction_pads*2.);
      hist_diff->GetXaxis()->SetTitleSize(hist_diff->GetXaxis()->GetTitleSize()/(fraction_pads));
      hist_diff->GetXaxis()->SetTitleOffset((hist_diff->GetXaxis()->GetTitleOffset()-(1.0-fraction_pads))/(fraction_pads));
      hist_diff->GetYaxis()->SetNdivisions(503);
      hist_diff->GetYaxis()->SetLabelSize(hist_diff->GetYaxis()->GetLabelSize()/(fraction_pads));
      hist_diff->GetYaxis()->SetLabelOffset(hist_diff->GetYaxis()->GetLabelOffset()/(fraction_pads));
      hist_diff->GetYaxis()->SetTitleSize(hist_diff->GetYaxis()->GetTitleSize()/(fraction_pads));
      hist_diff->GetYaxis()->SetTitleOffset(hist_diff->GetYaxis()->GetTitleOffset()*(fraction_pads));
      hist_diff->SetStats(kFALSE);

      graph_diff_exp = new TGraphAsymmErrors(nbins);
      graph_diff_exp->SetLineWidth(2);
      graph_diff_exp->SetMarkerStyle(0);
      graph_diff_exp->SetFillColor(kYellow);
      for (int i = 0; i < nbins; ++i) {
         hist_diff->SetBinContent(i+1, fHistData.GetBinContent(i+1)-histsum->GetBinContent(i+1));
         hist_diff->SetBinError(i+1, fHistData.GetBinError(i+1));
         graph_diff_exp->SetPoint(i, (graph_error_exp->GetX())[i], 0.0);
         graph_diff_exp->SetPointEXlow(i, 0.0);
         graph_diff_exp->SetPointEXhigh(i, 0.0);
         graph_diff_exp->SetPointEYlow(i, (graph_error_exp->GetEYlow())[i]);
         graph_diff_exp->SetPointEYhigh(i, (graph_error_exp->GetEYhigh())[i]);

         if (-(graph_error_exp->GetEYlow())[i] < ymin)
            ymin = -(graph_error_exp->GetEYlow())[i];
         if ((graph_error_exp->GetEYhigh())[i] > ymax)
            ymax = (graph_error_exp->GetEYhigh())[i];
      }
      if (ymax < (hist_diff->GetMaximum() + hist_diff->GetBinError(hist_diff->GetMaximumBin())))
         ymax = 1.1 * (hist_diff->GetMaximum() + hist_diff->GetBinError(hist_diff->GetMaximumBin()));
      if (ymin>(hist_diff->GetMinimum() - hist_diff->GetBinError(hist_diff->GetMaximumBin())))
         ymin = 1.1 * (hist_diff->GetMinimum() - hist_diff->GetBinError(hist_diff->GetMaximumBin()));
      (hist_diff->GetYaxis())->SetRangeUser(TMath::Floor(-1.1*TMath::Max(-ymin, ymax)), TMath::Ceil(1.1*TMath::Max(-ymin, ymax)));

   }

   // draw histograms
   stack.Draw("SAMEAF");
   stack.GetHistogram() -> SetXTitle("");
   stack.GetHistogram() -> SetYTitle("");
   stack.GetHistogram() -> GetXaxis() -> SetLabelSize(0);
   stack.GetHistogram() -> GetYaxis() -> SetLabelSize(0);
   stack.Draw("SAME");
   fHistData.Draw("SAMEP");

   if (flag_error0)
      fHistData.Draw("SAMEPE");

   if (flag_error1)
      histsum->Draw("SAMEE");

   if (flag_error3)
      graph_error_obs->Draw("SAMEZ");

   if (flag_error2)
      graph_error_exp->Draw("SAME||");

   if (flag_legend) {
      legend1->Draw();
      if (ntemplates>2)
      legend2->Draw();
   }

   TLine * line = 0;
   if (flag_diff) {
      pad2->cd();
      hist_diff->Draw("P");
      graph_diff_exp->Draw("SAME4");
      line = new TLine((hist_diff->GetXaxis())->GetXmin(), 0.0, (hist_diff->GetXaxis())->GetXmax(), 0.0);
      line->SetLineWidth(2);
      line->SetLineColor(kBlack);
      line->Draw("SAME");
      hist_diff->Draw("SAMEP");
   }

   c1->Print(filename);

   // delete temporary histograms
   delete pad1;
   delete pad2;
    delete c1;
   delete legend1;
    delete legend2;
    delete graph_error_exp;
    delete graph_error_obs;
    delete histsum;
   if (flag_diff) {
       delete hist_diff;
       delete graph_diff_exp;
       delete line;
   }
}

void BCTemplateFitter::PrintTemp

(

)

Temporary entry. Used for debugging.

Definition at line 1288 of file BCTemplateFitter.cxx.

{
   TCanvas * c1 = new TCanvas("");

   c1->cd();
   fUncertaintyHistogramExp->Draw("COL");
   c1->Print("uncertainty_exp.eps");

   c1->cd();
   fUncertaintyHistogramObsPosterior->Draw("COL");
   c1->Print("uncertainty_obs.eps");

   delete c1;
}

int BCTemplateFitter::PrintTemplate	(	const char *	name,
		const char *	filename
	)

Print a template to a file.

Parameters:

	name	The template name.
	filename	The filename.

Returns:

An error code.

Definition at line 1304 of file BCTemplateFitter.cxx.

{
   // get number of sources of systematic uncertainty
   int nsyst = GetNSystErrors();

   // get index
   int index = GetIndexTemplate(name);

   // check index
   if (index < 0) {
      BCLog::OutError("BCTemplateFitter::PrintTemplate : Could not find template.");
      return 0;
   }

   // remember old directory
   TDirectory* dir = gDirectory;

   // create postscript
   TPostScript * ps = new TPostScript(filename);

   // create new canvas
   TCanvas * c1 = new TCanvas("", "", 700, 700);

   c1->Update();
   ps->NewPage();
   c1->cd();

   // create legend
   TLegend l1(0.18, 0.75, 0.85, 0.85);
   l1.SetBorderSize(0);
   l1.SetFillColor(kWhite);

   // draw histogram and uncertainties
   TH1D hist_template = fTemplateHistogramContainer.at(index);
   hist_template.SetFillColor(kWhite);
   hist_template.SetFillStyle(0);
   hist_template.SetMarkerSize(0);
   hist_template.SetLineWidth(0);
   l1.AddEntry(&hist_template, name, "L");
   TH1D hist_totalerr = CombineUncertainties(name);
   TH1D hist_template_totalerr = CreateErrorHist(hist_template, hist_totalerr);
   hist_template_totalerr.SetFillColor(kYellow);
   hist_template_totalerr.SetFillStyle(1001);
   hist_template_totalerr.SetMarkerSize(0);
   l1.AddEntry(&hist_template_totalerr, "Systematic uncertainties", "F");
   TH1D hist_efferr = fEffErrHistogramContainer.at(index);
   TH1D hist_template_efferr = CreateErrorHist(hist_template, hist_efferr);
   hist_template_totalerr.SetFillColor(kRed);
   hist_template_totalerr.SetFillStyle(1001);
   hist_template_efferr.SetMarkerSize(0);
   l1.AddEntry(&hist_template_efferr, "Efficiency uncertainties", "F");
   int binmax = hist_template.GetMaximumBin();
   double ymax = hist_template.GetBinContent(binmax) + 2.0 * hist_template_totalerr.GetBinError(binmax);
   hist_template_totalerr.GetYaxis()->SetRangeUser(0.0, 1.25 * ymax);
   hist_template_totalerr.Draw("E2");
   hist_template_efferr.Draw("SAMEE2");
   hist_template.Draw("SAME");

   // draw legend
   l1.Draw();

   // update ps
   c1->Update();
   ps->NewPage();
   c1->cd();

   // create legend
   TLegend l2(0.18, 0.75, 0.85, 0.85);
   l2.SetBorderSize(0);
   l2.SetFillColor(kWhite);

   // print uncertainties
   c1->cd(2);
   hist_efferr = fEffErrHistogramContainer.at(index);
   double ymin = hist_efferr.GetMinimum();
   ymax = hist_efferr.GetMaximum();
   l2.AddEntry(&hist_efferr, "Efficiency", "L");
   hist_efferr.SetStats(kFALSE);
   hist_efferr.Draw();

   // loop over all uncertainties
   for (int i = 0; i < nsyst; ++i) {
      TH1D * hist = new TH1D(fSystErrorHistogramContainer.at(i).at(index));
      hist->SetLineColor(2 + i);
      if (hist->GetMaximum()>ymax)
         ymax = hist->GetMaximum();
      if (hist->GetMinimum()<ymin)
         ymin = hist->GetMinimum();
      l2.AddEntry(hist, fSystErrorNameContainer.at(i).c_str(), "L");
      hist->Draw("SAME");
   }
   if (ymin < 0)
      ymin = 1.25*ymin;
   else
      ymin = 0.8*ymin;

   if (ymax > 0)
      ymax = 1.25*ymax;
   else
      ymax = 0.8*ymax;

   hist_efferr.GetYaxis()->SetRangeUser(ymin, ymax);

   // draw legend
   l2.Draw();

   // close ps
   c1->Update();
   ps->Close();

   // change to old directory
   dir->cd();

   // free memory
   delete c1;
   delete ps;

   // no error
   return 1;
}

int BCTemplateFitter::SetData

(

const TH1D &

hist

)

Set the histogram containing the data.

Parameters:

hist

The data histogram.

Returns:

An error code.

Definition at line 222 of file BCTemplateFitter.cxx.

{
   // compare previous data set with this one
   if (fHistData.Integral() > 0.) {
      if (CompareHistogramProperties(fHistData, hist) != 1) {
         BCLog::OutError("BCTemplateFitter::SetData : Data and this histogram properties are incompatible.");
         return 0;
      }
   }

   // create histogram
   fNBins = hist.GetNbinsX();
   fXmin = (hist.GetXaxis())->GetXmin();
   fXmax = (hist.GetXaxis())->GetXmax();
   fHistData = TH1D(hist);

   // copy histogram content
   for (int i = 1; i <= fNBins; ++i)
      fHistData.SetBinContent(i, hist.GetBinContent(i));

   // set histogram style
   fHistData.SetXTitle((hist.GetXaxis())->GetTitle());
   fHistData.SetYTitle((hist.GetYaxis())->GetTitle());
   fHistData.SetMarkerStyle(20);
   fHistData.SetMarkerSize(1.1);
   fHistData.SetStats(kFALSE);

   // calculate norm
   fNorm = hist.Integral();

   // no errors
   return 1;
}

void BCTemplateFitter::SetFlagFixNorm

(

bool

flag

)

[inline]

Set the flag for using a fixed normalization (true) or floating normalization (false).

Definition at line 223 of file BCTemplateFitter.h.

         { fFlagFixNorm = flag; }

void BCTemplateFitter::SetFlagPhysicalLimits

(

bool

flag

)

[inline]

Set a flag for having physical limits (expectation values greater or equal to 0).

Parameters:

flag

The flag.

Definition at line 216 of file BCTemplateFitter.h.

         { fFlagPhysicalLimits = flag; }

void BCTemplateFitter::SetNorm

(

double

norm

)

[inline]

Set normalization constant.

Definition at line 229 of file BCTemplateFitter.h.

         { fNorm = norm; }

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

Describe the efficiency and the uncertainty for all bins.

Parameters:

	name	The template name.
	effmean	The mean value of the efficiency.
	errsigma	The uncertainty on the efficiency.

Definition at line 1247 of file BCTemplateFitter.cxx.

{
   // get index
   int index = GetIndexTemplate(name);

   // check index
   if (index < 0) {
      BCLog::OutError("BCTemplateFitter::SetTemplateEfficiency : Could not find template.");
      return 0;
   }

   // create histograms
   TH1D histeff = TH1D(fHistData);
    TH1D histefferr = TH1D(fHistData);

   // fill histograms
   for (int i = 1; i <= fNBins; ++i) {
      histeff.SetBinContent(i, effmean);
      histefferr.SetBinContent(i, effsigma);
   }

   // set histograms
   int err = SetTemplateEfficiency(name, histeff, histefferr);

   // return error code
   return err;
}

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

Describe the efficiency and the uncertainty for each bin.

Parameters:

	name	The template name.
	eff	A histogram describing the efficieny.
	efferr	A histogram describing the uncertainty on the efficiency.

Returns:

An error code.

Definition at line 1213 of file BCTemplateFitter.cxx.

{
   // get index
   int index = GetIndexTemplate(name);

   // check index
   if (index < 0) {
      BCLog::OutError("BCTemplateFitter::SetTemplateEfficiency : Could not find template.");
      return 0;
   }

   // check efficiency histogram
   if (fTemplateTypeContainer[index]==0)
      if (CompareHistogramProperties(fTemplateHistogramContainer.at(index), eff) != 1) {
         BCLog::OutError("BCTemplateFitter::SetTemplate efficiency : Template and efficiency histogram properties are incompatible.");
         return 0;
      }

   // set histogram style
   eff.SetXTitle((fHistData.GetXaxis())->GetTitle());
   eff.SetYTitle("Efficiency");

   efferr.SetXTitle((fHistData.GetXaxis())->GetTitle());
   efferr.SetYTitle("Efficiency uncertainty");

   // set efficiency histogram
   fEffHistogramContainer[index] = eff;
   fEffErrHistogramContainer[index] = efferr;

   // no error
   return 1;
}

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

Set the systematic uncertainty on a template.

Parameters:

	errorname	The name of the source.
	templatename	The template name.
	parerror	A histogram describing the uncertainty.

Returns:

An error code.

Definition at line 459 of file BCTemplateFitter.cxx.

{
   // get error index
   int errindex = GetIndexSystError(errorname);

   // check parameter index
   if (errindex < 0) {
      BCLog::OutError("BCTemplateFitter::SetTemplateSystError : Did not find parameter.");
      return 0;
   }

   // get template index
   int tempindex = GetIndexTemplate(templatename);

   // check index
   if (tempindex < 0) {
      BCLog::OutError("BCTemplateFitter::SetTemplateSystError : Could not find template.");
      return 0;
   }

   // set style
   parerror.SetStats(kFALSE);

   // set histogram
   (fSystErrorHistogramContainer.at(errindex))[tempindex] = parerror;

   // no error
   return 1;
}

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

Return the efficiency of a template in a certain bin.

Parameters:

	templatenumber	The number of the template.
	binnumber	The bin number.
	parameters	The parameter values.

Definition at line 134 of file BCTemplateFitter.cxx.

{
   // get number of sources of systematic uncertainties
   int nsyst = GetNSystErrors();

   // get efficiency index
   int effindex = fEffParIndexContainer.at(templatenumber);

   // get efficiency for the bin
   double efficiency = fEffHistogramContainer.at(templatenumber).GetBinContent(binnumber);

   // modify efficiency by uncertainty
   double efferr = fEffErrHistogramContainer.at(templatenumber).GetBinContent(binnumber);

   // check that efficiency error is larger than 0
   // add efficiency contribution from template
   efficiency += parameters.at(effindex) * efferr;

   // loop over sources of systematic uncertainties
   for (int isyst = 0; isyst < nsyst; ++isyst) {
      // get parameter index
      int systindex = fSystErrorParIndexContainer.at(isyst);

      // add efficiency
      double deff = fSystErrorHistogramContainer.at(isyst).at(templatenumber).GetBinContent(binnumber);

      // add efficiency contribution from systematic
      efficiency += parameters.at(systindex) * deff;
   }

   // make sure efficiency is between 0 and 1
   if (efficiency < 0.)
      efficiency = 0.;
   if (efficiency > 1.)
      efficiency = 1.;

   // return efficiency
   return efficiency;
}

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

Return the probability of a template in a certain bin.

Parameters:

	templateindex	The number of the template.
	binnumber	The bin number.
	parameters	The parameter values.

Definition at line 175 of file BCTemplateFitter.cxx.

{
   // initialize probability
   double probability = 0;

   // get template type
   int templatetype = fTemplateTypeContainer[templatenumber];

   // calculate probability from template histogram
   if (templatetype == 0) {
      int index = GetHistogramIndex(templatenumber);

      probability = fTemplateHistogramContainer.at(index).GetBinContent(binnumber);
   }
   else if (templatetype == 1) {
      double bincenter = fHistData.GetBinCenter(binnumber);
      int index = GetFunctionIndex(templatenumber);
      probability = fTemplateFunctionContainer.at(index).Eval(bincenter);
   }

   // return probability
   return probability;
}

---

Member Data Documentation

std::vector< std::vector<int> > BCTemplateFitter::fConstraintSumIndices [protected]

A container for constrained sums: indices.

Definition at line 502 of file BCTemplateFitter.h.

std::vector< double > BCTemplateFitter::fConstraintSumMean [protected]

A container for constrained sums: mean values.

Definition at line 507 of file BCTemplateFitter.h.

std::vector< double > BCTemplateFitter::fConstraintSumRMS [protected]

A container for constrained sums: mean values.

Definition at line 512 of file BCTemplateFitter.h.

std::vector<TH1D> BCTemplateFitter::fEffErrHistogramContainer [protected]

A container of efficiency uncertainty histograms.

Definition at line 447 of file BCTemplateFitter.h.

std::vector<TH1D> BCTemplateFitter::fEffHistogramContainer [protected]

A container of efficiency histograms.

Definition at line 443 of file BCTemplateFitter.h.

std::vector<int> BCTemplateFitter::fEffParIndexContainer [protected]

A container of parameter indeces for efficiencies.

Definition at line 479 of file BCTemplateFitter.h.

bool BCTemplateFitter::fFlagFixNorm [protected]

Flag for fixing the normalization or not

Definition at line 532 of file BCTemplateFitter.h.

bool BCTemplateFitter::fFlagPhysicalLimits [protected]

Flag for having physical limits (expectation values greater or equal to 0).

Definition at line 538 of file BCTemplateFitter.h.

std::vector<int> BCTemplateFitter::fFunctionIndexContainer [protected]

A container of function indeces for templates.

Definition at line 475 of file BCTemplateFitter.h.

TH1D BCTemplateFitter::fHistData [protected]

The data histogram.

Definition at line 429 of file BCTemplateFitter.h.

TH1D BCTemplateFitter::fHistNorm [protected]

Histogram containing the overall number of expected events.

Definition at line 517 of file BCTemplateFitter.h.

std::vector<int> BCTemplateFitter::fHistogramIndexContainer [protected]

A container of histogram indeces for templates.

Definition at line 471 of file BCTemplateFitter.h.

std::vector<TH1D> BCTemplateFitter::fHistRatios1D [protected]

1-D histograms containing ratios.

Definition at line 527 of file BCTemplateFitter.h.

std::vector< std::vector<int> > BCTemplateFitter::fIndicesRatios1D [protected]

Vector of indices for the calculation of ratios.

Definition at line 522 of file BCTemplateFitter.h.

int BCTemplateFitter::fNBins [protected]

The number of bins in the data.

Definition at line 557 of file BCTemplateFitter.h.

double BCTemplateFitter::fNorm [protected]

Normalization constant

Definition at line 553 of file BCTemplateFitter.h.

std::vector< std::vector<TH1D> > BCTemplateFitter::fSystErrorHistogramContainer [protected]

A matrix of histograms describing systematic uncertainties.

Definition at line 451 of file BCTemplateFitter.h.

std::vector<std::string> BCTemplateFitter::fSystErrorNameContainer [protected]

A container of names of sources of systematic uncertainties.

Definition at line 461 of file BCTemplateFitter.h.

std::vector<int> BCTemplateFitter::fSystErrorParIndexContainer [protected]

A container of parameter indeces for sources of systematic uncertainty.

Definition at line 484 of file BCTemplateFitter.h.

std::vector<std::string> BCTemplateFitter::fSystErrorTypeContainer [protected]

A container of error types.

Definition at line 490 of file BCTemplateFitter.h.

std::vector<TF1> BCTemplateFitter::fTemplateFunctionContainer [protected]

A container of template functions.

Definition at line 439 of file BCTemplateFitter.h.

std::vector<TH1D> BCTemplateFitter::fTemplateHistogramContainer [protected]

A container of template histograms.

Definition at line 435 of file BCTemplateFitter.h.

std::vector<std::string> BCTemplateFitter::fTemplateNameContainer [protected]

A container of template names.

Definition at line 457 of file BCTemplateFitter.h.

std::vector<int> BCTemplateFitter::fTemplateParIndexContainer [protected]

A container of parameter indeces for templates.

Definition at line 467 of file BCTemplateFitter.h.

std::vector<int> BCTemplateFitter::fTemplateTypeContainer [protected]

A container for template types.
0: histogram 1: function

Definition at line 497 of file BCTemplateFitter.h.

TH2D* BCTemplateFitter::fUncertaintyHistogramExp [protected]

A 2-d histogram for calculating the error bars

Definition at line 543 of file BCTemplateFitter.h.

TH2D* BCTemplateFitter::fUncertaintyHistogramObsPosterior [protected]

A 2-d histogram for calculating the error bars

Definition at line 548 of file BCTemplateFitter.h.

double BCTemplateFitter::fXmax [protected]

The maximum value of the data range.

Definition at line 565 of file BCTemplateFitter.h.

double BCTemplateFitter::fXmin [protected]

The minimum value of the data range.

Definition at line 561 of file BCTemplateFitter.h.

---

The documentation for this class was generated from the following files:

• BCTemplateFitter.h
• BCTemplateFitter.cxx