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)
	BCTemplateFitter ()
	BCTemplateFitter (const char *name)
double	CalculateChi2 ()
double	CalculateChi2Prob ()
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)
TH1D	GetData ()
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	GetParIndexEff (const char *name)
int	GetParIndexSystError (const char *name)
int	GetParIndexTemplate (const char *name)
int	GetParIndexTemplate (int index)
int	Initialize ()
double	LogAPrioriProbability (std::vector< double > parameters)
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, TH1D eff, TH1D efferr)
int	SetTemplateEfficiency (const char *name, double effmean=1., double effsigma=0.)
int	SetTemplatePrior (const char *name, TH1D prior)
int	SetTemplatePrior (const char *name, double mean, double sigma)
int	SetTemplateSystError (const char *errorname, const char *templatename, TH1D parerror)
	~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
TH1D	fHistData
TH1D	fHistNorm
std::vector< TH1D >	fHistRatios1D
std::vector< std::vector< int > >	fIndicesRatios1D
int	fNBins
double	fNorm
std::vector< TH1D >	fPriorContainer
int	fPriorNBins
std::vector< std::vector< TH1D > >	fSystErrorHistogramContainer
std::vector< std::string >	fSystErrorNameContainer
std::vector< int >	fSystErrorParIndexContainer
std::vector< std::string >	fSystErrorTypeContainer
std::vector< TH1D >	fTemplateHistogramContainer
std::vector< std::string >	fTemplateNameContainer
std::vector< int >	fTemplateParIndexContainer
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 34 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.)
   , fPriorNBins(1000)
{
}

BCTemplateFitter::BCTemplateFitter

(

const char *

name

)

A constructor.

Definition at line 38 of file BCTemplateFitter.cxx.

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

BCTemplateFitter::~BCTemplateFitter

(

)

The default destructor.

Definition at line 53 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 368 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("systerr_%i", GetNSystErrors()), -dx, dx);

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

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

   // 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 218 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 (CompareHistogramProperties(fHistData, hist) != 1) {
      BCLog::OutError("BCTemplateFitter::AddTemplate : Data and template histogram properties are incompatible.");
      return 0;
   }

   // 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 = int(fTemplateHistogramContainer.size());

   // 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 histprior = TH1D("", "", fPriorNBins, Nmin, Nmax);
   TH1D histsysterror = TH1D("", "", fNBins, fXmin, fXmax);

   // set style
   histprior.SetXTitle(name);

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

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

   // add a parameter for the expectation value
   AddParameter(Form("N_%i", partemplateindex), Nmin, Nmax);

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

   // add a parameter for the efficiency
   AddParameter(Form("eff_%i", partemplateindex), -5.0, 5.0);

   // add histogram, name and index to containers
   fTemplateHistogramContainer.push_back(hist);
   fPriorContainer.push_back(histprior);
   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.);

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

double BCTemplateFitter::CalculateChi2

(

)

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

Definition at line 871 of file BCTemplateFitter.cxx.

{
   int nbins = fHistData.GetNbinsX();
   int ntemplates = int(fTemplateHistogramContainer.size());

   std::vector <double> parameters = GetBestFitParameters();

   double chi2 = 0;

   // loop over all bins
   for (int ibin = 1; ibin <= nbins; ++ibin) {
      double nexp = 0;
      double ndata = fHistData.GetBinContent(ibin);

      // 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 + parameters.at(effindex) * efferr);

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

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

   // return chi2
   return chi2;
}

double BCTemplateFitter::CalculateChi2Prob

(

)

Calculates and returns the chi2-probability.

Definition at line 913 of file BCTemplateFitter.cxx.

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

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

double BCTemplateFitter::CalculateKSProb

(

)

Calculates and returns the Kolmogorov-Smirnov-probability.

Definition at line 930 of file BCTemplateFitter.cxx.

{
   // create a histogram with the sum of all contributions
   TH1 * histsum = (TH1D*)(fTemplateHistogramContainer.at(0)).Clone("temp");

   int nbins = fHistData.GetNbinsX();
   int ntemplates = int(fTemplateHistogramContainer.size());

   std::vector <double> parameters = GetBestFitParameters();

   // loop over all bins
   for (int ibin = 1; ibin <= nbins; ++ibin) {
      double bincontent = 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 + parameters.at(effindex) * efferr);

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

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

   // 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 923 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 977 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;
   }

   // define temporary variables
   int nbins = fHistData.GetNbinsX();
   int ntemplates = int( fTemplateHistogramContainer.size() );

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

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

   // define starting distribution
   for (int ibin = 1; ibin <= nbins; ++ibin) {
      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 < nbins; ++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 499 of file BCTemplateFitter.cxx.

{
   // get number of templates
   int ntemplates = int( fTemplateHistogramContainer.size() );

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

   // 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 544 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 1228 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 1378 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;
}

TH1D BCTemplateFitter::GetData

(

)

[inline]

Return a template histogram.

Parameters:

name

The template name. Return the histogram containing the data.

Definition at line 134 of file BCTemplateFitter.h.

         { return fHistData; }

TH1D BCTemplateFitter::GetHistRatio1D

(

int

index

)

[inline]

Return a ratio histogram

Definition at line 86 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 80 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 1464 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 1445 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 68 of file BCTemplateFitter.h.

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

int BCTemplateFitter::GetNRatios

(

)

[inline]

Return the number of ratios to calculate.

Definition at line 74 of file BCTemplateFitter.h.

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

int BCTemplateFitter::GetNSystErrors

(

)

[inline]

Return the number of sources of systematic uncertainties.

Definition at line 62 of file BCTemplateFitter.h.

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

int BCTemplateFitter::GetNTemplates

(

)

[inline]

Return the number of templates.

Definition at line 56 of file BCTemplateFitter.h.

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

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 1517 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 1536 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 1483 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 1502 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);
}

int BCTemplateFitter::Initialize

(

)

Initialize the fitting procedure.

Returns:

An error code.

Definition at line 449 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();

   fUncertaintyHistogramExp = new TH2D(
         TString::Format("UncertaintyExp_%i", BCLog::GetHIndex()), "",
         fHistData.GetNbinsX(), fHistData.GetXaxis()->GetXmin(), fHistData.GetXaxis()->GetXmax(),
         100, 0., maximum);

   fUncertaintyHistogramObsPosterior = new TH2D(
         TString::Format("UncertaintyObsPosterior_%i", BCLog::GetHIndex()), "",
         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 = int(fTemplateParIndexContainer.size());

   // 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::LogAPrioriProbability

(

std::vector< double >

parameters

)

[virtual]

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

Reimplemented from BCModel.

Definition at line 132 of file BCTemplateFitter.cxx.

{
   double logprob = 0.;

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

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

      // prior on process contributions
      double par = parameters.at(fTemplateParIndexContainer.at(i));
      int bin = fPriorContainer.at(i).FindBin(par);
      logprob += log( fPriorContainer.at(i).GetBinContent(bin) );

      // prior on efficiences
      par = parameters.at(fEffParIndexContainer.at(i));
      logprob += BCMath::LogGaus(par, 0.0, 1.0);
   }

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

   // loop over sources of systematic uncertainties
   for (int i = 0; i < nsyst; ++i) {
      double par = parameters.at(fSystErrorParIndexContainer.at(i));
      if (fSystErrorTypeContainer.at(i) == "gauss")
         logprob += BCMath::LogGaus(par, 0.0, 1.0);
   }

   // constraints
   int nconstraints = int(fConstraintSumIndices.size());
   if (nconstraints > 0) {

      // loop over constraints
      for (int i = 0; i < nconstraints; ++i) {

         // initialize sum
         double sum = 0;

         // get number of summands
         int nsummands = int( (fConstraintSumIndices.at(i)).size() );

         // loop over summands and add to sum
         for (int j = 0; j < nsummands; ++j) {
            int index = fConstraintSumIndices.at(i).at(j);
            double par = parameters.at(fTemplateParIndexContainer.at(index));
            sum += par;
         }

         // add to prior
         logprob += BCMath::LogGaus(sum, fConstraintSumMean.at(i), fConstraintSumRMS.at(i));
      }
   }

   return logprob;
}

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 63 of file BCTemplateFitter.cxx.

{
   double logprob = 0.;

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

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

   // loop over bins
   for (int ibin = 1; ibin <= fNBins; ++ibin) {
      double nexp = 0;
      double ndata = fHistData.GetBinContent(ibin);

      // loop over all templates
      for (int itemp = 0; itemp < ntemplates; ++itemp) {
         int templateindex = 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 += 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(itemp).GetBinContent(ibin);

            if (deff > 0)
               efficiency += deff * parameters.at(systindex);
         }

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

         // calculate expectation nvalue
         nexp += parameters.at(templateindex)
               * efficiency
               * fTemplateHistogramContainer.at(itemp).GetBinContent(ibin);
      }

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

      // 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 1077 of file BCTemplateFitter.cxx.

{
   int nbins      = fHistData.GetNbinsX();
   int ntemplates = int(fTemplateHistogramContainer.size());

   // loop over all bins
   for (int ibin = 1; ibin <= nbins; ++ibin) {
      double bincontent = 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 + fMCMCx.at(effindex) * efferr);

         bincontent += fMCMCx.at(tempindex) * efficiency * fTemplateHistogramContainer.at(itemp).GetBinContent(ibin);
      }

      // 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 1555 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 1137 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();

   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 563 of file BCTemplateFitter.cxx.

{
   int nbins = fHistData.GetNbinsX();
   int ntemplates = int( fTemplateHistogramContainer.size() );

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

   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;

   // 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);
      pad1->Draw();
      pad2->Draw();
   }
   else {
      pad1 = new TPad("pad1", "",0.0, 0.0, 1.0, 1.0);
      pad1->SetFillColor(kWhite);
      pad2 = new TPad();
      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");

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

   // scale histograms and add to stack and legend
   for (int itemp = 0; itemp < ntemplates; ++itemp) {
      int tempindex = fTemplateParIndexContainer.at(itemp);
      int effindex = fEffParIndexContainer.at(itemp);

      // scale histogram
      fTemplateHistogramContainer.at(itemp).Scale(
            GetBestFitParameter(tempindex)/ fTemplateHistogramContainer.at(itemp).Integral());

      // loop over bins and scale these
      for (int ibin = 1; ibin <= fNBins; ++ibin) {
         // 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 + GetBestFitParameter(effindex) * efferr);

         fTemplateHistogramContainer.at(itemp).SetBinContent(
               ibin,fTemplateHistogramContainer.at(itemp).GetBinContent(ibin) * efficiency);
      }

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

   // loop over all bins
   for (int ibin = 1; ibin <= nbins; ++ibin) {
      double bincontent = 0;

      // loop over all templates
      for (int itemp = 0; itemp < ntemplates; ++itemp)
         bincontent +=fTemplateHistogramContainer.at(itemp).GetBinContent(ibin);

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

   // 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->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(-1.1*TMath::Max(-ymin, ymax), 1.1*TMath::Max(-ymin, ymax));

   }

   // draw histograms
   stack.Draw("SAMEA");
   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);

   // rescale
   for (int i = 0; i < ntemplates; ++i)
      fTemplateHistogramContainer.at(i).Scale(1.0 / fTemplateHistogramContainer.at(i).Integral());

   // 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 1240 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 1256 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 191 of file BCTemplateFitter.cxx.

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

   // 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 149 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 142 of file BCTemplateFitter.h.

         { fFlagPhysicalLimits = flag; }

void BCTemplateFitter::SetNorm

(

double

norm

)

[inline]

Set normalization constant.

Definition at line 155 of file BCTemplateFitter.h.

         { fNorm = norm; }

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 1166 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 (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::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 1199 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("", "", fNBins, fXmin, fXmax);
   TH1D histefferr = TH1D("", "", fNBins, fXmin, fXmax);

   // 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::SetTemplatePrior	(	const char *	name,
		double	mean,
		double	sigma
	)

Set a Gaussian prior on the template.

Parameters:

	name	The name of the template.
	mean	The mean value of the prior.
	rms	The standard deviation of the prior.

Returns:

An error code.

Definition at line 336 of file BCTemplateFitter.cxx.

{
   // get index
   int parindex = GetParIndexTemplate(name);

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

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

   // create new histogram
   TH1D hist("", "", fPriorNBins, par->GetLowerLimit(), par->GetUpperLimit());

   // loop over bins and fill histogram
   for (int i = 1; i < fPriorNBins; ++i) {
      double x = hist.GetBinCenter(i);
      double fx = TMath::Gaus(x, mean, sigma);
      hist.SetBinContent(i, fx);
   }

   // set template prior
   int err = SetTemplatePrior(name, hist);

   // return error code
   return err;
}

int BCTemplateFitter::SetTemplatePrior	(	const char *	name,
		TH1D	prior
	)

Set an arbitrary prior on the template

Parameters:

	name	The name of the template.
	prior	A histogram describing the prior.

Returns:

An error code.

Definition at line 308 of file BCTemplateFitter.cxx.

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

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

   // check prior
   if (prior.Integral() <= 0) {
      BCLog::OutError("BCTemplateFitter::SetTemplatePrior : Integral of prior is equal to zero or less than that.");
      return 0;
   }

   // normalize prior to unity
   prior.Scale(1.0/prior.Integral());

   // replace histogram
   fPriorContainer[parindex] = prior;

   // 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 418 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;
}

---

Member Data Documentation

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

A container for constrained sums: indices.

Definition at line 432 of file BCTemplateFitter.h.

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

A container for constrained sums: mean values.

Definition at line 437 of file BCTemplateFitter.h.

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

A container for constrained sums: mean values.

Definition at line 442 of file BCTemplateFitter.h.

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

A container of efficiency uncertainty histograms.

Definition at line 388 of file BCTemplateFitter.h.

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

A container of efficiency histograms.

Definition at line 384 of file BCTemplateFitter.h.

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

A container of parameter indeces for efficiencies.

Definition at line 416 of file BCTemplateFitter.h.

bool BCTemplateFitter::fFlagFixNorm [protected]

Flag for fixing the normalization or not

Definition at line 462 of file BCTemplateFitter.h.

bool BCTemplateFitter::fFlagPhysicalLimits [protected]

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

Definition at line 468 of file BCTemplateFitter.h.

TH1D BCTemplateFitter::fHistData [protected]

The data histogram.

Definition at line 374 of file BCTemplateFitter.h.

TH1D BCTemplateFitter::fHistNorm [protected]

Histogram containing the overall number of expected events.

Definition at line 447 of file BCTemplateFitter.h.

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

1-D histograms containing ratios.

Definition at line 457 of file BCTemplateFitter.h.

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

Vector of indices for the calculation of ratios.

Definition at line 452 of file BCTemplateFitter.h.

int BCTemplateFitter::fNBins [protected]

The number of bins in the data.

Definition at line 487 of file BCTemplateFitter.h.

double BCTemplateFitter::fNorm [protected]

Normalization constant

Definition at line 483 of file BCTemplateFitter.h.

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

A container of prior histograms.

Reimplemented from BCModel.

Definition at line 396 of file BCTemplateFitter.h.

int BCTemplateFitter::fPriorNBins [protected]

The number of bins for a prior distribution.

Definition at line 499 of file BCTemplateFitter.h.

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

A matrix of histograms describing systematic uncertainties.

Definition at line 392 of file BCTemplateFitter.h.

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

A container of names of sources of systematic uncertainties.

Definition at line 406 of file BCTemplateFitter.h.

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

A container of parameter indeces for sources of systematic uncertainty.

Definition at line 421 of file BCTemplateFitter.h.

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

A container of error types.

Definition at line 427 of file BCTemplateFitter.h.

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

A container of template histograms.

Definition at line 380 of file BCTemplateFitter.h.

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

A container of template names.

Definition at line 402 of file BCTemplateFitter.h.

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

A container of parameter indeces for templates.

Definition at line 412 of file BCTemplateFitter.h.

TH2D* BCTemplateFitter::fUncertaintyHistogramExp [protected]

A 2-d histogram for calculating the error bars

Definition at line 473 of file BCTemplateFitter.h.

TH2D* BCTemplateFitter::fUncertaintyHistogramObsPosterior [protected]

A 2-d histogram for calculating the error bars

Definition at line 478 of file BCTemplateFitter.h.

double BCTemplateFitter::fXmax [protected]

The maximum value of the data range.

Definition at line 495 of file BCTemplateFitter.h.

double BCTemplateFitter::fXmin [protected]

The minimum value of the data range.

Definition at line 491 of file BCTemplateFitter.h.

---

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

• BCTemplateFitter.h
• BCTemplateFitter.cxx