BCTemplateFitter Class Reference

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

#include <BCTemplateFitter.h>

Inheritance diagram for BCTemplateFitter:

Collaboration diagram for BCTemplateFitter:

List of all members.

Public Member Functions
int	AddSystError (const char *errorname, const char *errtype="gauss")
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="")
	BCTemplateFitter (const char *name)
	BCTemplateFitter ()
double	CalculateChi2 (std::vector< double > parameters)
double	CalculateChi2Prob (std::vector< double > parameters)
double	CalculateKSProb ()
double	CalculateMaxLike ()
double	CalculatePValue ()
int	CalculateRatio (int index, std::vector< int > indices, double rmin=-1.0, double rmax=1.0)
TH1D	CombineUncertainties (const char *name)
int	CompareHistogramProperties (TH1D hist1, TH1D hist2)
int	ConstrainSum (std::vector< int > indices, double mean, double rms)
TH1D	CreateErrorHist (TH1D hist, TH1D histerr)
double	Expectation (int binnumber, std::vector< double > &parameters)
int	FixTemplateFunctions (std::vector< double > &parameters)
TH1D	GetData ()
int	GetFunctionIndex (int templatenumber)
int	GetHistogramIndex (int templatenumber)
TH1D	GetHistRatio1D (int index)
std::vector< TH1D >	GetHistRatios1D ()
int	GetIndexSystError (const char *name)
int	GetIndexTemplate (const char *name)
int	GetNDF ()
int	GetNRatios ()
int	GetNSystErrors ()
int	GetNTemplates ()
int	GetNTemplatesType (int templatetype)
int	GetParIndexEff (const char *name)
int	GetParIndexSystError (const char *name)
int	GetParIndexTemplate (int index)
int	GetParIndexTemplate (const char *name)
std::vector< TF1 >	GetTemplateFunctionContainer ()
std::vector< TH1D >	GetTemplateHistogramContainer ()
std::vector< int >	GetTemplateParIndexContainer ()
int	Initialize ()
double	LogLikelihood (std::vector< double > parameters)
void	MCMCUserIterationInterface ()
int	PerformFit ()
void	PrintRatios (const char *filename="ratio.ps", int option=0, double ovalue=0.)
void	PrintStack (const char *filename="stack.ps", const char *options="LE2E0D")
void	PrintTemp ()
int	PrintTemplate (const char *name, const char *filename)
int	SetData (const TH1D &hist)
void	SetFlagFixNorm (bool flag)
void	SetFlagPhysicalLimits (bool flag)
void	SetNorm (double norm)
int	SetTemplateEfficiency (const char *name, double effmean=1., double effsigma=0.)
int	SetTemplateEfficiency (const char *name, TH1D eff, TH1D efferr)
int	SetTemplateSystError (const char *errorname, const char *templatename, TH1D parerror)
double	TemplateEfficiency (int templatenumber, int binnumber, std::vector< double > &parameters)
double	TemplateProbability (int templatenumber, int binnumber, std::vector< double > &parameters)
	~BCTemplateFitter ()
Protected Attributes
std::vector< std::vector< int > >	fConstraintSumIndices
std::vector< double >	fConstraintSumMean
std::vector< double >	fConstraintSumRMS
std::vector< TH1D >	fEffErrHistogramContainer
std::vector< TH1D >	fEffHistogramContainer
std::vector< int >	fEffParIndexContainer
bool	fFlagFixNorm
bool	fFlagPhysicalLimits
std::vector< int >	fFunctionIndexContainer
TH1D	fHistData
TH1D	fHistNorm
std::vector< int >	fHistogramIndexContainer
std::vector< TH1D >	fHistRatios1D
std::vector< std::vector< int > >	fIndicesRatios1D
int	fNBins
double	fNorm
std::vector< std::vector< TH1D > >	fSystErrorHistogramContainer
std::vector< std::string >	fSystErrorNameContainer
std::vector< int >	fSystErrorParIndexContainer
std::vector< std::string >	fSystErrorTypeContainer
std::vector< TF1 >	fTemplateFunctionContainer
std::vector< TH1D >	fTemplateHistogramContainer
std::vector< std::string >	fTemplateNameContainer
std::vector< int >	fTemplateParIndexContainer
std::vector< int >	fTemplateTypeContainer
TH2D *	fUncertaintyHistogramExp
TH2D *	fUncertaintyHistogramObsPosterior
double	fXmax
double	fXmin

---

Detailed Description

A class for fitting several templates to a data set.

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 37 of file BCTemplateFitter.h.

---

Constructor & Destructor Documentation

BCTemplateFitter::BCTemplateFitter

(

)

The default constructor.

Definition at line 23 of file BCTemplateFitter.cxx.

   : BCModel()
   , fFlagFixNorm(false)
   , fFlagPhysicalLimits(true)
   , fUncertaintyHistogramExp(0)
   , fUncertaintyHistogramObsPosterior(0)
   , fNorm(-1)
   , fNBins(-1)
   , fXmin(1.)
   , fXmax(0.)
{
}

BCTemplateFitter::BCTemplateFitter

(

const char *

name

)

A constructor.

Definition at line 37 of file BCTemplateFitter.cxx.

   : BCModel(name)
   , fFlagFixNorm(false)
   , fFlagPhysicalLimits(true)
   , fUncertaintyHistogramExp(0)
   , fUncertaintyHistogramObsPosterior(0)
   , fNorm(-1)
   , fNBins(-1)
   , fXmin(1.)
   , fXmax(0.)
{
}

BCTemplateFitter::~BCTemplateFitter

(

)

The default destructor.

Definition at line 51 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 396 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 247 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 340 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

(

std::vector< double >

parameters

)

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

Definition at line 971 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

(

std::vector< double >

parameters

)

Calculates and returns the chi2-probability.

Definition at line 992 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 1009 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 1002 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 1036 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 540 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 1450 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 585 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 1272 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 1422 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,
		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 85 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

(

std::vector< double > &

parameters

)

Definition at line 190 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 155 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 169 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 162 of file BCTemplateFitter.h.

      { return fHistogramIndexContainer.at(templatenumber); }; 

TH1D BCTemplateFitter::GetHistRatio1D

(

int

index

)

[inline]

Return a ratio histogram

Definition at line 107 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 89 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 1509 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 1490 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 77 of file BCTemplateFitter.h.

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

int BCTemplateFitter::GetNRatios

(

)

[inline]

Return the number of ratios to calculate.

Definition at line 83 of file BCTemplateFitter.h.

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

int BCTemplateFitter::GetNSystErrors

(

)

[inline]

Return the number of sources of systematic uncertainties.

Definition at line 71 of file BCTemplateFitter.h.

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

int BCTemplateFitter::GetNTemplates

(

)

[inline]

Return the number of templates.

Definition at line 59 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 1619 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 1562 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 1581 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 1528 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 1547 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 101 of file BCTemplateFitter.h.

         { return fTemplateFunctionContainer; }; 

std::vector<TH1D> BCTemplateFitter::GetTemplateHistogramContainer

(

)

[inline]

Return the container of template histograms.

Definition at line 95 of file BCTemplateFitter.h.

         { return fTemplateHistogramContainer; }; 

std::vector<int> BCTemplateFitter::GetTemplateParIndexContainer

(

)

[inline]

Return container of parameter indeces for templates

Definition at line 175 of file BCTemplateFitter.h.

         { return fTemplateParIndexContainer; };

int BCTemplateFitter::Initialize

(

)

Initialize the fitting procedure.

Returns:

An error code.

Definition at line 480 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());
   
   Double_t 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,
         //       fHistData.GetNbinsX(), fHistData.GetXaxis()->GetXmin(), fHistData.GetXaxis()->GetXmax(),
         100, minimum, maximum);

   fUncertaintyHistogramObsPosterior = new TH2D(
         TString::Format("UncertaintyObsPosterior_%i", BCLog::GetHIndex()), "",
         fHistData.GetNbinsX(), a,
         //       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

(

std::vector< double >

parameters

)

[virtual]

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

Reimplemented from BCModel.

Definition at line 61 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 1139 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 1600 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 1179 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 604 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);
      Double_t 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);
      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 1284 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 1300 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 212 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 214 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 207 of file BCTemplateFitter.h.

         { fFlagPhysicalLimits = flag; }

void BCTemplateFitter::SetNorm

(

double

norm

)

[inline]

Set normalization constant.

Definition at line 220 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 1243 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 1209 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 449 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,
		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 121 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
   if (efferr > 0)
      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);
      
      // check that efficiency error is larger than 0
      // add efficiency contribution from systematic
      if (deff > 0)
         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,
		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 165 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 493 of file BCTemplateFitter.h.

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

A container for constrained sums: mean values.

Definition at line 498 of file BCTemplateFitter.h.

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

A container for constrained sums: mean values.

Definition at line 503 of file BCTemplateFitter.h.

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

A container of efficiency uncertainty histograms.

Definition at line 438 of file BCTemplateFitter.h.

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

A container of efficiency histograms.

Definition at line 434 of file BCTemplateFitter.h.

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

A container of parameter indeces for efficiencies.

Definition at line 470 of file BCTemplateFitter.h.

bool BCTemplateFitter::fFlagFixNorm [protected]

Flag for fixing the normalization or not

Definition at line 523 of file BCTemplateFitter.h.

bool BCTemplateFitter::fFlagPhysicalLimits [protected]

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

Definition at line 529 of file BCTemplateFitter.h.

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

A container of function indeces for templates.

Definition at line 466 of file BCTemplateFitter.h.

TH1D BCTemplateFitter::fHistData [protected]

The data histogram.

Definition at line 420 of file BCTemplateFitter.h.

TH1D BCTemplateFitter::fHistNorm [protected]

Histogram containing the overall number of expected events.

Definition at line 508 of file BCTemplateFitter.h.

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

A container of histogram indeces for templates.

Definition at line 462 of file BCTemplateFitter.h.

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

1-D histograms containing ratios.

Definition at line 518 of file BCTemplateFitter.h.

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

Vector of indices for the calculation of ratios.

Definition at line 513 of file BCTemplateFitter.h.

int BCTemplateFitter::fNBins [protected]

The number of bins in the data.

Definition at line 548 of file BCTemplateFitter.h.

double BCTemplateFitter::fNorm [protected]

Normalization constant

Definition at line 544 of file BCTemplateFitter.h.

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

A matrix of histograms describing systematic uncertainties.

Definition at line 442 of file BCTemplateFitter.h.

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

A container of names of sources of systematic uncertainties.

Definition at line 452 of file BCTemplateFitter.h.

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

A container of parameter indeces for sources of systematic uncertainty.

Definition at line 475 of file BCTemplateFitter.h.

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

A container of error types.

Definition at line 481 of file BCTemplateFitter.h.

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

A container of template functions.

Definition at line 430 of file BCTemplateFitter.h.

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

A container of template histograms.

Definition at line 426 of file BCTemplateFitter.h.

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

A container of template names.

Definition at line 448 of file BCTemplateFitter.h.

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

A container of parameter indeces for templates.

Definition at line 458 of file BCTemplateFitter.h.

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

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

Definition at line 488 of file BCTemplateFitter.h.

TH2D* BCTemplateFitter::fUncertaintyHistogramExp [protected]

A 2-d histogram for calculating the error bars

Definition at line 534 of file BCTemplateFitter.h.

TH2D* BCTemplateFitter::fUncertaintyHistogramObsPosterior [protected]

A 2-d histogram for calculating the error bars

Definition at line 539 of file BCTemplateFitter.h.

double BCTemplateFitter::fXmax [protected]

The maximum value of the data range.

Definition at line 556 of file BCTemplateFitter.h.

double BCTemplateFitter::fXmin [protected]

The minimum value of the data range.

Definition at line 552 of file BCTemplateFitter.h.

---

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

• BCTemplateFitter.h
• BCTemplateFitter.cxx