Private Attributes

BCMTF Class Reference

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

#include <BCMTF.h>

Inherits BCModel.

Collaboration diagram for BCMTF:
Collaboration graph
[legend]

List of all members.

Public Member Functions

Constructors and destructors

 BCMTF ()
 BCMTF (const char *name)
 ~BCMTF ()
Member functions (get)

int GetNChannels ()
int GetNProcesses ()
int GetNSystematics ()
int GetChannelIndex (const char *name)
int GetProcessIndex (const char *name)
int GetSystematicIndex (const char *name)
int GetParIndexProcess (int index)
int GetParIndexSystematic (int index)
BCMTFChannel * GetChannel (int index)
BCMTFProcess * GetProcess (int index)
BCMTFSystematic * GetSystematic (int index)
Member functions (set)

int SetData (const char *channelname, TH1D hist, double minimum=-1, double maximum=-1)
int SetTemplate (const char *channelname, const char *processname, TH1D hist, double efficiency=1.)
int SetTemplate (const char *channelname, const char *processname, std::vector< TF1 * > *funccont, int nbins, double efficiency=1.)
void SetExpectationFunction (int parindex, TF1 *func)
int SetSystematicVariation (const char *channelname, const char *processname, const char *systematicname, double variation_up, double variation_down)
int SetSystematicVariation (const char *channelname, const char *processname, const char *systematicname, TH1D hist_up, TH1D hist_down)
int SetSystematicVariation (const char *channelname, const char *processname, const char *systematicname, TH1D hist, TH1D hist_up, TH1D hist_down)
void SetFlagEfficiencyConstraint (bool flag)
Member functions (miscellaneous methods)

int AddChannel (const char *name)
int AddProcess (const char *name, double nmin=0., double nmax=1.)
int AddSystematic (const char *name, double min=-5., double max=5.)
double Expectation (int channelindex, int binindex, const std::vector< double > &parameters)
double ExpectationFunction (int parindex, int channelindex, int processindex, const std::vector< double > &parameters)
double Efficiency (int channelindex, int processindex, int binindex, const std::vector< double > &parameters)
double Probability (int channelindex, int processindex, int binindex, const std::vector< double > &parameters)
double CalculateChi2 (int channelindex, const std::vector< double > &parameters)
double CalculateChi2 (const std::vector< double > &parameters)
double CalculateCash (int channelindex, const std::vector< double > &parameters)
double CalculateCash (const std::vector< double > &parameters)
Member functions (output methods)

int PrintSummary (const char *filename="summary.txt")
int PrintStack (int channelindex, const std::vector< double > &parameters, const char *filename="stack.eps", const char *options="e1b0stack")
int PrintStack (const char *channelname, const std::vector< double > &parameters, const char *filename="stack.eps", const char *options="e1b0stack")
Member functions (overloaded from BCModel)

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

Private Attributes

std::vector< BCMTFChannel * > fChannelContainer
std::vector< BCMTFProcess * > fProcessContainer
std::vector< BCMTFSystematic * > fSystematicContainer
int fNChannels
int fNProcesses
int fNSystematics
std::vector< int > fProcessParIndexContainer
std::vector< int > fSystematicParIndexContainer
bool fFlagEfficiencyConstraint
std::vector< TF1 * > fExpectationFunctionContainer

Detailed Description

A class for fitting several templates to a data set.

Author:
Daniel Kollar
Kevin Kröninger
Version:
1.1
Date:
06.2012 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.

Definition at line 37 of file BCMTF.h.

Constructor & Destructor Documentation

BCMTF::BCMTF (  ) 

The default constructor.

Definition at line 34 of file BCMTF.cxx.

 : BCModel("Multi-template Fitter")
 , fNChannels(0)
 , fNProcesses(0)
 , fNSystematics(0)
 , fFlagEfficiencyConstraint(false)
{}

BCMTF::BCMTF ( const char *  name  ) 

A constructor.

Parameters:
name The name of the model

Definition at line 43 of file BCMTF.cxx.

 : BCModel(name)
 , fNChannels(0)
 , fNProcesses(0)
 , fNSystematics(0)
 , fFlagEfficiencyConstraint(false)
{}

BCMTF::~BCMTF (  ) 

The default destructor.

Definition at line 52 of file BCMTF.cxx.

{
   for (int i = 0; i < fNChannels; ++i)
      delete fChannelContainer.at(i);
}

Member Function Documentation

int BCMTF::AddChannel ( const char *  name  ) 

Add a channel

Parameters:
name The channel name.
Returns:
An error code.

Definition at line 274 of file BCMTF.cxx.

{
   // check if channel exists
   for (int i = 0; i < fNChannels; ++i) {
      // compare names
      if (GetChannelIndex(name) >= 0) {
      BCLog::OutWarning("BCMultitemplateFitter::AddChannel() : Channel with this name exists already.");
         return -1;
      }
   }

   // create new channel
   BCMTFChannel * channel = new BCMTFChannel(name);

   // create new data
   BCMTFTemplate * bctemplate = new BCMTFTemplate(channel->GetName().c_str(), "data");

   // add data
   channel->SetData(bctemplate);

   // add process templates
   for (int i = 0; i < fNProcesses; ++i) {
      // get process
      BCMTFProcess * process = GetProcess(i);

      // create new template
      BCMTFTemplate * bctemplate = new BCMTFTemplate(name, process->GetName().c_str());

      // add template
      channel->AddTemplate(bctemplate);
   }

   // loop over all systematics
   for (int i = 0; i < fNSystematics; ++i) {
      // get systematic
      BCMTFSystematic * systematic = GetSystematic(i);

      // create new systematic variation
      BCMTFSystematicVariation * variation = new BCMTFSystematicVariation(name, systematic->GetName().c_str(), fNProcesses);

      // add systematic variation
      channel->AddSystematicVariation(variation);
   }

   // add channel
   fChannelContainer.push_back(channel);

   // increase number of channels
   fNChannels++;

   // no error
   return 1;
}

int BCMTF::AddProcess ( const char *  name,
double  nmin = 0.,
double  nmax = 1. 
)

Add a process and the associated BAT parameter.

Parameters:
name The process name
nmin The minimum number of expected events (lower limit on the BAT parameter values).
nmax The maximum number of expected events (upper limit on the BAT parameter values).
Returns:
An error code.

Definition at line 329 of file BCMTF.cxx.

{
   // check if process exists
   for (int i = 0; i < fNProcesses; ++i) {
      // compare names
      if (GetProcessIndex(name) >= 0) {
      BCLog::OutWarning("BCMultitemplateFitter::AddProcess() : Process with this name exists already.");
         return -1;
      }
   }

   // create new process
   BCMTFProcess * process = new BCMTFProcess(name);

   // add process
   fProcessContainer.push_back(process);

   // add process templates
   for (int i = 0; i < fNChannels; ++i) {
      // get channel
      BCMTFChannel * channel = GetChannel(i);

      // create new template
      BCMTFTemplate * bctemplate = new BCMTFTemplate(channel->GetName().c_str(), name);

      // add template
      channel->AddTemplate(bctemplate);

      // loop over all systematic
      for (int j = 0; j < fNSystematics; ++j) {
         // get systematic variation
         BCMTFSystematicVariation * variation = channel->GetSystematicVariation(j);

         // add histogram
         variation->AddHistograms(0, 0);
      }
   }

   // increase number of processes
   fNProcesses++;

   // add parameter index to container
   fProcessParIndexContainer.push_back(GetNParameters());

   // add parameter
   AddParameter(name, nmin, nmax);

   // add a functional form for the expectation
   fExpectationFunctionContainer.push_back(0);

   // no error
   return 1;
}

int BCMTF::AddSystematic ( const char *  name,
double  min = -5.,
double  max = 5. 
)

Add a source of systematic uncertainty and the associated BAT (nuisance) parameter.

Parameters:
name The systematic uncertainty name.
min The lower limit on the BAT parameter values, typically -5 sigma if Gaussian constraint is used.
max The upper limit on the BAT parameter values, typically +5 sigma if Gaussian constraint is used.
Returns:

Definition at line 384 of file BCMTF.cxx.

{
   // check if systematic exists
   for (int i = 0; i < fNSystematics; ++i) {
      // compare names
      if (GetSystematicIndex(name) >= 0) {
      BCLog::OutWarning("BCMultitemplateFitter::AddSystematic() : Systematic with this name exists already.");
         return -1;
      }
   }

   // create new systematic
   BCMTFSystematic * systematic = new BCMTFSystematic(name);

   // add systematic
   fSystematicContainer.push_back(systematic);

   // add systematic variations
   for (int i = 0; i < fNChannels; ++i) {
      // get channel
      BCMTFChannel * channel = GetChannel(i);

      // create new systematic variation
      BCMTFSystematicVariation * variation = new BCMTFSystematicVariation(channel->GetName().c_str(), name, fNProcesses);

      // add systematic variation
      channel->AddSystematicVariation(variation);
   }
   // ...

   // increase number of systematices
   fNSystematics++;

   // add parameter index to container
   fSystematicParIndexContainer.push_back(GetNParameters());

   // add parameter
   AddParameter(name, min, max);

   // add a functional form for the expectation
   fExpectationFunctionContainer.push_back(0);

   // no error
   return 1;
}

double BCMTF::CalculateCash ( int  channelindex,
const std::vector< double > &  parameters 
)

Calculate the Cash statistic for a single channel

Parameters:
channelindex The channel index.
parameters A reference to the parameters used to calculate the Cash statistic.
Returns:
The Cash statistic.
See also:
CalculateCash(const std::vector<double> & parameters)

Definition at line 1049 of file BCMTF.cxx.

{
   if (parameters.size() == 0)
      return -1;

   double cash = 0;

   // get channel
   BCMTFChannel * channel = GetChannel(channelindex);

   // get data
   BCMTFTemplate * data = channel->GetData();

   // get histogram
   TH1D * hist = data->GetHistogram();

   // check if histogram exists
   if (hist) {
      // get number of bins in data
      int nbins = hist->GetNbinsX();

      // loop over all bins
      for (int ibin = 1; ibin <= nbins; ++ibin) {
         // get expectation value
         double expectation = Expectation(channelindex, ibin, parameters);

         // get observation
         double observation = hist->GetBinContent(ibin);

         // calculate Cash statistic
         cash += 2. * (expectation - observation);

         // check negative log
         if (observation > 0)
            cash += 2. * observation * log (observation/expectation);
      }
   }

   // return cash;
   return cash;

}

double BCMTF::CalculateCash ( const std::vector< double > &  parameters  ) 

Calculate the Cash statistic for all channels

Parameters:
parameters A reference to the parameters used to calculate the Cash statistic.
Returns:
The Cash statistic.
See also:
CalculateCash(int channelindex, const std::vector<double> & parameters)

Definition at line 1093 of file BCMTF.cxx.

{
   if (parameters.size() == 0)
      return -1;

   double cash = 0;

   // get number of channels
   int nchannels = GetNChannels();

   // loop over all channels
   for (int i = 0; i < nchannels; ++i) {
      cash += CalculateCash(i, parameters);
   }

   // return cash
   return cash;
}

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

Calculate a chi2 for all channels together given a set of parameters.

Parameters:
parameters A reference to the parameters used to calculate the chi2.
Returns:
A chi2 value.
See also:
CalculateChi2(int channelindex, const std::vector<double> & parameters)

Definition at line 1029 of file BCMTF.cxx.

{
   if (parameters.size() == 0)
      return -1;

   double chi2 = 0;

   // get number of channels
   int nchannels = GetNChannels();

   // loop over all channels
   for (int i = 0; i < nchannels; ++i) {
      chi2 += CalculateChi2(i, parameters);
   }

   // return chi2
   return chi2;
}

double BCMTF::CalculateChi2 ( int  channelindex,
const std::vector< double > &  parameters 
)

Calculate a chi2 for a single channel given a set of parameters.

Parameters:
channelindex The channel index.
parameters A reference to the parameters used to calculate the chi2.
Returns:
A chi2 value.
See also:
CalculateChi2(const std::vector<double> & parameters)

Definition at line 990 of file BCMTF.cxx.

{
   if (parameters.size() == 0)
      return -1;

   double chi2 = 0;

   // get channel
   BCMTFChannel * channel = GetChannel(channelindex);

   // get data
   BCMTFTemplate * data = channel->GetData();

   // get histogram
   TH1D * hist = data->GetHistogram();

   // check if histogram exists
   if (hist) {
      // get number of bins in data
      int nbins = hist->GetNbinsX();

      // loop over all bins
      for (int ibin = 1; ibin <= nbins; ++ibin) {
         // get expectation value
         double expectation = Expectation(channelindex, ibin, parameters);

         // get observation
         double observation = hist->GetBinContent(ibin);

         // add Poisson term
         chi2 += (expectation - observation) * (expectation - observation) / expectation;
      }
   }

   // return chi2;
   return chi2;
}

double BCMTF::Efficiency ( int  channelindex,
int  processindex,
int  binindex,
const std::vector< double > &  parameters 
)

Return the efficiency for a process in a channel and for a particular bin.

Parameters:
channelindex The channel index.
processindex The process index.
binindex The bin index.
parameters A reference to the parameters used to calculate the efficiency.
Returns:
The efficiency.

Definition at line 678 of file BCMTF.cxx.

{
   // get channel
   BCMTFChannel * channel = fChannelContainer[channelindex];

   double efficiency = channel->GetTemplate(processindex)->GetEfficiency();

   double defficiency = 1.;

   // loop over all systematics
   for (int i = 0; i < fNSystematics; ++i) {
      if (!(fSystematicContainer[i]->GetFlagSystematicActive()))
         continue;

      // get parameter index
      int parindex = fSystematicParIndexContainer[i];

      // get histogram
      TH1D * hist = 0;

      if (parameters[parindex] > 0)
         hist = channel->GetSystematicVariation(i)->GetHistogramUp(processindex);
      else
         hist = channel->GetSystematicVariation(i)->GetHistogramDown(processindex);

      // check if histogram exists
      if (!hist)
         continue;

      // multiply efficiency
      defficiency += parameters[parindex] * hist->GetBinContent(binindex);
   }

   // calculate efficiency
   efficiency *= defficiency;

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

   return efficiency;
}

double BCMTF::Expectation ( int  channelindex,
int  binindex,
const std::vector< double > &  parameters 
)

Return the expected number of events for a channel and bin.

Parameters:
channelindex The channel index.
binindex The bin index.
parameters A reference to the parameters used to calculate the expectation.
Returns:
The expectation value

Definition at line 631 of file BCMTF.cxx.

{
   double expectation = 0.;
    
   // loop over all processes
   for (int i = 0; i < fNProcesses; ++i) {
       // get efficiency
       double efficiency = Efficiency(channelindex, i, binindex, parameters);
       
       // get probability
       double probability = Probability(channelindex, i, binindex, parameters);
       
       // get parameter index
       int parindex = fProcessParIndexContainer[i];
       
       // add to expectation
       expectation += ExpectationFunction(parindex, channelindex, i, parameters)
          * efficiency
          * probability;
   }
    
   // check if expectation is positive
   if (expectation < 0)
      expectation = 0.;

   return expectation;
}

double BCMTF::ExpectationFunction ( int  parindex,
int  channelindex,
int  processindex,
const std::vector< double > &  parameters 
)

Return the function value of the expectation function for a parameter, channel and process.

Parameters:
parindex The parameter index.
channelindex The channel index.
processindex The process index.
A reference to the parameters used to calculate the expectation.
Returns:
The expectation function value.

Definition at line 660 of file BCMTF.cxx.

{
  // get function container
  std::vector<TF1 *> * funccont = fChannelContainer[channelindex]->GetTemplate(processindex)->GetFunctionContainer();

  if (funccont->size()>0)
    return 1.;

  else if (!fExpectationFunctionContainer[parindex])
    return parameters[parindex];

  else {
    TF1 * func = fExpectationFunctionContainer[parindex];
    return func->Eval(parameters[parindex]);
  }
}

BCMTFChannel* BCMTF::GetChannel ( int  index  )  [inline]
Parameters:
name The channel index.
Returns:
The channel object.

Definition at line 107 of file BCMTF.h.

         { return fChannelContainer.at(index); };

int BCMTF::GetChannelIndex ( const char *  name  ) 
Parameters:
name The name of the channel.
Returns:
The channel index.

Definition at line 60 of file BCMTF.cxx.

{
   // loop over all channels and compare names
   for (int i = 0; i < fNChannels; ++i) {
      // get channel
      BCMTFChannel * channel = GetChannel(i);

      // compare names
      if (!channel->GetName().compare(name))
         return i;
   }

   // if channel does not exist, return -1
   return -1;
}

int BCMTF::GetNChannels (  )  [inline]
Returns:
The number of channels.

Definition at line 64 of file BCMTF.h.

         { return fNChannels; };

int BCMTF::GetNProcesses (  )  [inline]
Returns:
The number of processes.

Definition at line 69 of file BCMTF.h.

         { return fNProcesses; };

int BCMTF::GetNSystematics (  )  [inline]
Returns:
The number of systematics.

Definition at line 74 of file BCMTF.h.

         { return fNSystematics; };

int BCMTF::GetParIndexProcess ( int  index  )  [inline]
Parameters:
index The parameter index (mtf counting) .
Returns:
The parameter number corresponding to the parameter index (BAT counting).

Definition at line 95 of file BCMTF.h.

         { return fProcessParIndexContainer.at(index); };

int BCMTF::GetParIndexSystematic ( int  index  )  [inline]
Parameters:
name The systematic uncertainty index (mtf counting).
Returns:
The parameter number corresponding to the systematic uncertainty (BAT counting).

Definition at line 101 of file BCMTF.h.

         { return fSystematicParIndexContainer.at(index); };

BCMTFProcess* BCMTF::GetProcess ( int  index  )  [inline]
Parameters:
name The process index.
Returns:
The process object.

Definition at line 113 of file BCMTF.h.

         { return fProcessContainer.at(index); };

int BCMTF::GetProcessIndex ( const char *  name  ) 
Parameters:
name The name of the process.
Returns:
The process index.

Definition at line 77 of file BCMTF.cxx.

{
   // loop over all processs and compare names
   for (int i = 0; i < fNProcesses; ++i) {
      // get process
      BCMTFProcess * process = GetProcess(i);

      // compare names
      if (!process->GetName().compare(name))
         return i;
   }

   // if process does not exist, return -1
   return -1;
}

BCMTFSystematic* BCMTF::GetSystematic ( int  index  )  [inline]
Parameters:
name The systematic ucnertainty index.
Returns:
The systematic uncertainty object.

Definition at line 119 of file BCMTF.h.

         { return fSystematicContainer.at(index); };

int BCMTF::GetSystematicIndex ( const char *  name  ) 
Parameters:
name The name of the systematic.
Returns:
The systematic uncertainty index.

Definition at line 94 of file BCMTF.cxx.

{
   // loop over all systematics and compare names
   for (int i = 0; i < fNSystematics; ++i) {
      // get systematic
      BCMTFSystematic * systematic = GetSystematic(i);

      // compare names
      if (!systematic->GetName().compare(name))
         return i;
   }

   // if process does not exist, return -1
   return -1;
}

double BCMTF::LogLikelihood ( const std::vector< double > &  parameters  )  [virtual]

Calculates natural logarithm of the likelihood. Method needs to be overloaded by the user.

Parameters:
params A set of parameter values
Returns:
Natural logarithm of the likelihood

Implements BCModel.

Definition at line 1113 of file BCMTF.cxx.

{
   double logprob = 0.;

   // loop over all channels
   for (int ichannel = 0; ichannel < fNChannels; ++ichannel) {

      // get channel
      BCMTFChannel * channel = fChannelContainer[ichannel];

      // check if channel is active
      if (!(channel->GetFlagChannelActive()))
         continue;

     // get data
      BCMTFTemplate * data = channel->GetData();

      // get histogram
      TH1D * hist = data->GetHistogram();

      // check if histogram exists
      if (!hist)
         continue;

      // get number of bins in data
      int nbins = data->GetNBins();

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

         // get expectation value
         double expectation = Expectation(ichannel, ibin, parameters);

         // get observation
         double observation = hist->GetBinContent(ibin);

         // add Poisson term
         logprob += BCMath::LogPoisson(observation, expectation);
      }
   }

   return logprob;
}

void BCMTF::MCMCUserIterationInterface (  )  [virtual]

Method executed for every iteration of the MCMC. User's code should be provided via overloading in the derived class

Reimplemented from BCIntegrate.

Definition at line 1158 of file BCMTF.cxx.

{
   // loop over all channels
   for (int ichannel = 0; ichannel < fNChannels; ++ichannel) {

      // get channel
      BCMTFChannel * channel = fChannelContainer[ichannel];

      // check if channel is active
      if (!(channel->GetFlagChannelActive()))
         continue;

     // get data
      BCMTFTemplate * data = channel->GetData();

      // get histogram
      TH1D * hist_data = data->GetHistogram();

      // check if histogram exists
      if (!hist_data)
         continue;

         // get histogram for uncertainty band
         TH2D* hist_uncbandexp = channel->GetHistUncertaintyBandExpectation(); 

      // check if histogram exists
      if (!hist_uncbandexp)
         continue;

      // get number of bins in data
      int nbins = hist_data->GetNbinsX();

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

         // get expectation value
         double expectation = Expectation(ichannel, ibin, fMCMCx);

             // fill uncertainty band on expectation
             hist_uncbandexp->Fill(hist_data->GetBinCenter(ibin), expectation);
      }
   }

}

int BCMTF::PrintStack ( int  channelindex,
const std::vector< double > &  parameters,
const char *  filename = "stack.eps",
const char *  options = "e1b0stack" 
)

Print the stack of templates together with the data in a particular channel. Several plot options are available:
"logx" : plot the x-axis on a log scale
"logy" : plot the y-axis on a log scale
"bw" : plot in black and white
"sum" : draw a line corresponding to the sum of all templates
"stack" : draw the templates as a stack
"e0" : do not draw error bars
"e1" : draw error bars corresponding to sqrt(n)
"b0" : draw an error band on the expectation corresponding to the central 68% probability
"b1" : draw bands showing the probability to observe a certain number of events given the expectation. The green (yellow, red) bands correspond to the central 68% (95%, 99.8%) probability

Parameters:
channelindex The channel index.
parameters A reference to the parameters used to scale the templates.
options The plotting options.
Returns:
An error code.

Definition at line 755 of file BCMTF.cxx.

{
   // todo:
   // - remove x-error on data points
   // - use hatched fill for error band

   // check if parameters are filled
   if (!parameters.size())
      return -1;

   // check options
   bool flag_logx   = false; // plot x-axis in log-scale
   bool flag_logy   = false; // plot y-axis in log-scale
   bool flag_bw     = false; // plot in black and white

    bool flag_sum    = false; // plot sum of all templates
    bool flag_stack  = false; // plot stack of templates

    bool flag_e0     = false; // do not draw error bars on data
    bool flag_e1     = false; // draw sqrt(N) error bars on data

    bool flag_b0     = false; // draw an error band on the expectation
    bool flag_b1     = false; // draw an error band on the number of events

   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;

   if (std::string(options).find("bw") < std::string(options).size())
      flag_bw = true;

   if (std::string(options).find("sum") < std::string(options).size())
      flag_sum = true;

   if (std::string(options).find("stack") < std::string(options).size())
      flag_stack = true;

   if (std::string(options).find("e0") < std::string(options).size())
      flag_e0 = true;

   if (std::string(options).find("e1") < std::string(options).size())
      flag_e1 = true;

   if (std::string(options).find("b0") < std::string(options).size())
      flag_b0 = true;

   if (std::string(options).find("b1") < std::string(options).size())
      flag_b1 = true;

    if (!flag_e0)
       flag_e1=true;

   // get channel
   BCMTFChannel * channel = GetChannel(channelindex);

   // create canvas
   TCanvas * c1 = new TCanvas();
   c1->cd();

   // set log or linear scale
   if (flag_logx)
      c1->SetLogx();

   if (flag_logy)
      c1->SetLogy();
    
    // get data histogram
    TH1D* hist_data = channel->GetData()->GetHistogram();

    // get number of bins
    int nbins = hist_data->GetNbinsX();

    // define sum of templates
    TH1D* hist_sum = new TH1D(*hist_data);
    hist_sum->SetLineColor(kBlack);
    for (int i = 1; i <= nbins; ++i)
       hist_sum->SetBinContent(i, 0);

  // define error band
    TH1D* hist_error_band = new TH1D(*hist_data);
    hist_error_band->SetFillColor(kBlack);
    hist_error_band->SetFillStyle(3005);
    hist_error_band->SetLineWidth(1);
    hist_error_band->SetStats(kFALSE);
    hist_error_band->SetMarkerSize(0);

   TGraphAsymmErrors * graph_error_exp = new TGraphAsymmErrors(nbins);
    //   graph_error_exp->SetLineWidth(2);
    graph_error_exp->SetMarkerStyle(0);
    graph_error_exp->SetFillColor(kBlack);
    graph_error_exp->SetFillStyle(3005);

    // get histogram for uncertainty band
    TH2D* hist_uncbandexp = channel->GetHistUncertaintyBandExpectation(); 

    // fill error band
  if (flag_b0) {
       for (int i = 1; i <= nbins; ++i) {
          TH1D * proj = hist_uncbandexp->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, hist_data->GetBinCenter(i), quantiles[1]);
          graph_error_exp->SetPointError(i-1, 0.0, 0.0, quantiles[1] - quantiles[0], quantiles[2]-quantiles[1]);
          hist_error_band->SetBinContent(i, 0.5*(quantiles[2]+quantiles[0]));
          hist_error_band->SetBinError(i, 0, 0.5*(quantiles[2]-quantiles[0]));
          delete proj;
       }
    }
    
   // create stack
   THStack * stack = new THStack("", "");

   // create a container of temporary histograms
   std::vector<TH1D *> histcontainer;

   // get number of templates
   unsigned int ntemplates = GetNProcesses();

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

      // get histogram
      TH1D * temphist = channel->GetTemplate(i)->GetHistogram();

      // get function container
      std::vector<TF1 *> * funccont = channel->GetTemplate(i)->GetFunctionContainer();

      // create new histogram
      TH1D * hist(0);

      if (temphist)
         hist = new TH1D( *(temphist) );
      else if (funccont)
         hist = new TH1D( *(channel->GetData()->GetHistogram()));
      else
         continue;

      if (flag_bw) {
         hist->SetFillColor(0);
         hist->SetLineWidth(1);
         hist->SetLineStyle(int(1+i));
      }
      else {
         hist->SetFillColor(2 + i);
         hist->SetFillStyle(1001);
      }

      // scale histogram
      for (int ibin = 1; ibin <= nbins; ++ibin) {

         // get efficiency
         double efficiency = Efficiency(channelindex, i, ibin, parameters);

         // get probability
         double probability = Probability(channelindex, i, ibin, parameters);

         // get parameter index
         int parindex = GetParIndexProcess(i);

         // add to expectation
         double expectation = parameters[parindex] * efficiency * probability;

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

             // add bin content
             hist_sum->SetBinContent(ibin, hist_sum->GetBinContent(ibin) + expectation);
      }

      // add histogram to container (for memory management)
      histcontainer.push_back(hist);

      // add histogram to stack
      stack->Add(hist);
   }

   //draw data
   hist_data->Draw("P0");

   // set range user
   hist_data->GetYaxis()->SetRangeUser(channel->GetRangeYMin(), channel->GetRangeYMax());

    // draw stack
    if (flag_stack)
       stack->Draw("SAMEHIST");

    // draw error band on number of observed events
    if (flag_b1) {
       channel->CalculateUncertaintyBandPoisson(0.001, 0.999, kRed)->Draw("SAMEE2");
       channel->CalculateUncertaintyBandPoisson(0.023, 0.977, kOrange)->Draw("SAMEE2");
       channel->CalculateUncertaintyBandPoisson(0.159, 0.841, kGreen)->Draw("SAMEE2");
    }

    // draw error band on expectation
    if (flag_b0) {
       hist_error_band->Draw("SAMEE2");
    }

    if (flag_sum)
       hist_sum->Draw("SAME");

    //draw data again
    if (flag_e0)
       hist_data->Draw("SAMEP0");

    if (flag_e1)
       hist_data->Draw("SAMEP0E");

   // redraw the axes
   gPad->RedrawAxis();

   // print
   c1->Print(filename);

   // free memory
   for (unsigned int i = 0; i < histcontainer.size(); ++i) {
      TH1D * hist = histcontainer.at(i);
      delete hist;
   }
   delete stack;
   delete c1;
    delete graph_error_exp;
    delete hist_error_band;
    delete hist_sum;

   // no error
   return 1;
}

int BCMTF::PrintStack ( const char *  channelname,
const std::vector< double > &  parameters,
const char *  filename = "stack.eps",
const char *  options = "e1b0stack" 
)

Print the stack of templates together with the data in a particular channel.

Parameters:
channelname The name of the channel.
parameters A reference to the parameters used to scale the templates.
options The plotting options.
Returns:
An error code.
See also:
PrintStack(int channelindex, const std::vector<double> & parameters, const char * filename = "stack.eps", const char * options = "e1b0stack")

Definition at line 747 of file BCMTF.cxx.

{
   int index = GetChannelIndex(channelname);

   return PrintStack(index, parameters, filename, options);
}

int BCMTF::PrintSummary ( const char *  filename = "summary.txt"  ) 

Print a summary of the fit into an ASCII file.

Parameters:
filename The name of the file.
Returns:
An error code

Definition at line 552 of file BCMTF.cxx.

{
   // open file
   std::ofstream ofi(filename);
   ofi.precision(3);

   // check if file is open
   if(!ofi.is_open()) {
      BCLog::OutWarning(Form("BCMultitemplateFitter::PrintSummary() : Could not open file %s", filename));
      return 0;
   }

   ofi
      << " Multi template fitter summary " << std::endl
      << " ----------------------------- " << std::endl
      << std::endl
      << " Number of channels      : " << fNChannels << std::endl
      << " Number of processes     : " << fNProcesses << std::endl
      << " Number of systematics   : " << fNSystematics << std::endl
      << std::endl;

   ofi
      << " Channels :" << std::endl;
      for (int i = 0; i < GetNChannels(); ++i) {
         ofi
            << " " << i
            << " : \"" << GetChannel(i)->GetName().c_str() << "\""
            << std::endl;
      }
      ofi
         << std::endl;

   ofi
      << " Processes :" << std::endl;
      for (int i = 0; i < GetNProcesses(); ++i) {
         ofi
            << " " << i
            << " : \"" << GetProcess(i)->GetName().c_str()  << "\""
            << " (par index " << GetParIndexProcess(i) << ")"
            << std::endl;
      }
      ofi
         << std::endl;

   ofi
      << " Systematics :" << std::endl;
      for (int i = 0; i < GetNSystematics(); ++i) {
         ofi
            << " " << i
            << " : \"" << GetSystematic(i)->GetName().c_str()  << "\""
            << " (par index " << GetParIndexSystematic(i) << ")"
            << std::endl;
      }
      ofi
         << std::endl;
      if (GetNSystematics() == 0)
         ofi
            << " - none - " << std::endl;

      ofi
         << " Goodness-of-fit: " << std::endl;
      for (int i = 0; i < GetNChannels(); ++i) {
         ofi
            << " i : \"" << GetChannel(i)->GetName().c_str() << "\" : chi2 = "
            << CalculateChi2( i, GetBestFitParameters() )
            << std::endl;
      }
      ofi
         << std::endl;


   // close file
   ofi.close();

   // no error
   return 1;
}

double BCMTF::Probability ( int  channelindex,
int  processindex,
int  binindex,
const std::vector< double > &  parameters 
)

Return the probability for a process in a channel and for a particular bin. This corresponds to the (normalized) bin content of the template.

Parameters:
channelindex The channel index.
processindex The process index.
binindex The bin index.
parameters A reference to the parameters used to calculate the probability.
Returns:
The probability.

Definition at line 726 of file BCMTF.cxx.

{
   // get histogram
   TH1D * hist = fChannelContainer[channelindex]->GetTemplate(processindex)->GetHistogram();

   // get function container
   std::vector<TF1 *> * funccont = fChannelContainer[channelindex]->GetTemplate(processindex)->GetFunctionContainer();

   // this needs to be fast
   if (!hist && !(funccont->size()>0))
      return 0.;

   if (hist)
     return hist->GetBinContent(binindex);
   else {
     int parindex = fProcessParIndexContainer[processindex];
     return funccont->at(binindex-1)->Eval(parameters[parindex]);
   }
}

int BCMTF::SetData ( const char *  channelname,
TH1D  hist,
double  minimum = -1,
double  maximum = -1 
)

Set the data histogram in a particular channel.

Parameters:
channelname The name of the channel.
hist The TH1D histogram.
minimum The minimum number of expected events (used for calculation of uncertainty bands).
maximum The maximum number of expected events (used for calculation of uncertainty bands).
Returns:
An error code.

Definition at line 199 of file BCMTF.cxx.

{
   int channelindex = GetChannelIndex(channelname);

   // check if channel exists
   if (channelindex < 0) {
      BCLog::OutWarning("BCMultitemplateFitter::SetTemplate() : Channel does not exist.");
         return -1;
   }

   // get channel
   BCMTFChannel * channel = GetChannel(channelindex);

   // get template
   BCMTFTemplate * data = channel->GetData();

   // remove statistics box
   hist.SetStats(kFALSE);

   // set marker
   hist.SetMarkerStyle(20);
    hist.SetMarkerSize(1.1);

    // set divisions
    hist.SetNdivisions(509);

   // remove old data set if it exists
   if (data->GetHistogram()) {
      delete data->GetHistogram();
      data->SetHistogram(0);
   }

    // remove old uncertainty histograms if they exist
    if (channel->GetHistUncertaintyBandExpectation()) {
       delete channel->GetHistUncertaintyBandExpectation();
       channel->SetHistUncertaintyBandExpectation(0);
    }
    if (channel->GetHistUncertaintyBandPoisson()) {
       delete channel->GetHistUncertaintyBandPoisson();
       channel->SetHistUncertaintyBandPoisson(0);
    }

    // create new histograms for uncertainty bands
    //    double minimum = floor(TMath::Max(0., hist.GetMinimum() - 7.*sqrt(hist.GetMinimum())));
    if (minimum==-1)
       minimum = 0;
    if (maximum==-1)
       maximum = ceil(hist.GetMaximum() + 5.*sqrt(hist.GetMaximum()));

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

   TH2D* hist_uncbandexp = new TH2D(TString::Format("UncertaintyBandExpectation_%i", BCLog::GetHIndex()), "",
                                                      hist.GetNbinsX(), &a[0], 1000, minimum, maximum);
    hist_uncbandexp->SetStats(kFALSE);

   TH2D* hist_uncbandpoisson = new TH2D(TString::Format("UncertaintyBandPoisson_%i", BCLog::GetHIndex()), "",
                                                            hist.GetNbinsX(), &a[0], int(maximum-minimum), minimum, maximum);
    hist_uncbandpoisson->SetStats(kFALSE);

   // set histograms
   data->SetHistogram(new TH1D(hist));
    channel->SetHistUncertaintyBandExpectation(hist_uncbandexp); 
    channel->SetHistUncertaintyBandPoisson(hist_uncbandpoisson); 

    // set y-range for printing
    channel->SetRangeY(minimum, maximum);

   // no error
   return 1;
}

void BCMTF::SetExpectationFunction ( int  parindex,
TF1 *  func 
) [inline]

Set an expectation function.

Parameters:
parindex The index of the parameter
func The pointer to a TF1 function.
See also:
SetTemplate(const char * channelname, const char * processname, std::vector<TF1 *> * funccont, int nbins, double efficiency = 1.)

Definition at line 166 of file BCMTF.h.

         { fExpectationFunctionContainer[parindex] = func; };

void BCMTF::SetFlagEfficiencyConstraint ( bool  flag  )  [inline]

Set a flag for the efficiency: if true then the total efficiency including all systematic uncertainties has to be between 0 and 1 at all times. Larger (smaller) values are set to 1 (0) during the calculation.

Parameters:
flag The flag

Definition at line 221 of file BCMTF.h.

         { fFlagEfficiencyConstraint = flag; };

int BCMTF::SetSystematicVariation ( const char *  channelname,
const char *  processname,
const char *  systematicname,
TH1D  hist_up,
TH1D  hist_down 
)

Set the impact of a source of systematic uncertainty for a particular source of systematic uncertainty and process in a given channel. The variation depends on x and is given as a histogram.

Parameters:
channelname The name of the channel.
processname The name of the process.
systematicname The name of the source of systematic uncertainty.
hist_up The TH1D histogram defining the relative shift between the up-variation and the nominal template: (up-nom)/nom.
hist_down The TH1D histogram defining the relative shift between the down-variation and the nominal template: (nom-down)/nom.
Returns:
An error code.
See also:
SetSystematicVariation(const char * channelname, const char * processname, const char * systematicname, double variation_up, double variation_down)

Definition at line 490 of file BCMTF.cxx.

{
   // get channel index
   int channelindex = GetChannelIndex(channelname);

   // check if channel exists
   if (channelindex < 0) {
      BCLog::OutWarning("BCMultitemplateFitter::SetTemplate() : Channel does not exist.");
         return -1;
   }

   // get process index
   int processindex = GetProcessIndex(processname);

   // check if process exists
   if (processindex < 0) {
      BCLog::OutWarning("BCMultitemplateFitter::SetTemplate() : Process does not exist.");
         return -1;
   }

   // get systematic index
   int systematicindex = GetSystematicIndex(systematicname);

   // check if systematic exists
   if (systematicindex < 0) {
      BCLog::OutWarning("BCMultitemplateFitter::SetTemplate() : Systematic does not exist.");
         return -1;
   }

   // get channel
   BCMTFChannel * channel = GetChannel(channelindex);

   // get systematic variation
   BCMTFSystematicVariation * variation = channel->GetSystematicVariation(systematicindex);

   // set histogram
   variation->SetHistograms(processindex, new TH1D(hist_up), new TH1D(hist_down));

   // no error
   return 1;
}

int BCMTF::SetSystematicVariation ( const char *  channelname,
const char *  processname,
const char *  systematicname,
TH1D  hist,
TH1D  hist_up,
TH1D  hist_down 
)

Set the impact of a source of systematic uncertainty for a particular source of systematic uncertainty and process in a given channel. The variation depends on x. The histograms are the raw histograms after the shift and the variation wrt the nominal template will be calculated automatically.

Parameters:
channelname The name of the channel.
processname The name of the process.
systematicname The name of the source of systematic uncertainty.
hist The histogram with the nominal template
hist_up The TH1D histogram after up-scaling of the systematic uncertainty.
hist_down The TH1D histogram after down-scaling of the systematic uncertainty.
Returns:
An error code.
See also:
SetSystematicVariation(const char * channelname, const char * processname, const char * systematicname, double variation_up, double variation_down)

Definition at line 533 of file BCMTF.cxx.

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

   TH1D * hist_up_rel = new TH1D(hist);
   TH1D * hist_down_rel = new TH1D(hist);

   // loop over all bins
   for (int ibin = 1; ibin <= nbins; ++ibin) {
      hist_up_rel->SetBinContent(ibin, (hist_up.GetBinContent(ibin) - hist.GetBinContent(ibin)) / hist.GetBinContent(ibin));
      hist_down_rel->SetBinContent(ibin, (hist.GetBinContent(ibin) - hist_down.GetBinContent(ibin)) / hist.GetBinContent(ibin));
   }

   // set the systematic variation
   return SetSystematicVariation(channelname, processname, systematicname, *hist_up_rel, *hist_down_rel);
}

int BCMTF::SetSystematicVariation ( const char *  channelname,
const char *  processname,
const char *  systematicname,
double  variation_up,
double  variation_down 
)

Set the impact of a source of systematic uncertainty for a particular source of systematic uncertainty and process in a given channel. The impact is the relative deviation between the varied and the nominal template, i.e., if the variation is set to 0.05 then the efficiency is multiplied by (1+0.05*eps), where eps is the corresponding nuisance parameter.

Parameters:
channelname The name of the channel.
processname The name of the process.
systematicname The name of the source of systematic uncertainty.
variation_up The relative shift between the up-variation and the nominal template: (up-nom)/nom.
variation_down The relative shift between the down-variation and the nominal template: (nom-down)/nom.
Returns:
An error code.

Definition at line 431 of file BCMTF.cxx.

{

   // get channel index
   int channelindex = GetChannelIndex(channelname);

   // check if channel exists
   if (channelindex < 0) {
      BCLog::OutWarning("BCMultitemplateFitter::SetTemplate() : Channel does not exist.");
         return -1;
   }

   // get process index
   int processindex = GetProcessIndex(processname);

   // check if process exists
   if (processindex < 0) {
      BCLog::OutWarning("BCMultitemplateFitter::SetTemplate() : Process does not exist.");
         return -1;
   }

   // get systematic index
   int systematicindex = GetSystematicIndex(systematicname);

   // check if systematic exists
   if (systematicindex < 0) {
      BCLog::OutWarning("BCMultitemplateFitter::SetTemplate() : Systematic does not exist.");
         return -1;
   }

   // get channel
   BCMTFChannel * channel = GetChannel(channelindex);

   BCMTFTemplate * bctemplate = channel->GetTemplate(processindex);

   TH1D * hist_template = bctemplate->GetHistogram();

   TH1D hist_up = TH1D(*hist_template);
   TH1D hist_down = TH1D(*hist_template);

   int nbins = hist_up.GetNbinsX();

   // loop over all bins
   for (int ibin = 1; ibin <= nbins; ++ibin) {
      hist_up.SetBinContent(ibin, variation_up);
      hist_down.SetBinContent(ibin, variation_down);
   }

   // get systematic variation
   BCMTFSystematicVariation * variation = channel->GetSystematicVariation(systematicindex);

   // set histogram
   variation->SetHistograms(processindex, new TH1D(hist_up), new TH1D(hist_down));

   // no error
   return 1;
}

int BCMTF::SetTemplate ( const char *  channelname,
const char *  processname,
TH1D  hist,
double  efficiency = 1. 
)

Set the template for a specific process in a particular channel.

Parameters:
channelname The name of the channel.
processname The name of the process.
hist The TH1D histogram.
efficiency The efficiency of this process in this channel.
Returns:
An error code.

Definition at line 111 of file BCMTF.cxx.

{
   // get channel index
   int channelindex = GetChannelIndex(channelname);

   // check if channel exists
   if (channelindex < 0) {
      BCLog::OutWarning("BCMultitemplateFitter::SetTemplate() : Channel does not exist.");
         return -1;
   }

   // get process index
   int processindex = GetProcessIndex(processname);

   // check if process exists
   if (processindex < 0) {
      BCLog::OutWarning("BCMultitemplateFitter::SetTemplate() : Process does not exist.");
         return -1;
   }

   // get channel
   BCMTFChannel * channel = GetChannel(channelindex);

   // get template
   BCMTFTemplate * bctemplate = channel->GetTemplate(processindex);

   // normalize histogram
   if (hist.Integral())
      hist.Scale(1.0 / hist.Integral());

   // remove statistics box
   hist.SetStats(kFALSE);

   // set color and fill style
   hist.SetFillColor(2 + processindex);
   hist.SetFillStyle(1001);

   // create new histogram
   TH1D * temphist = new TH1D(hist);

   // set histogram
   bctemplate->SetHistogram(temphist);

   // set efficiency
   bctemplate->SetEfficiency(efficiency);

   // no error
   return 1;
}

int BCMTF::SetTemplate ( const char *  channelname,
const char *  processname,
std::vector< TF1 * > *  funccont,
int  nbins,
double  efficiency = 1. 
)

Set the template for a specific process in a particular channel. This is an alternative way to describe the number of expected events. It is used in the rare case that processes cannot be summed directly, but interference effects have to be taken into account. The expected number of events is then parametrized as a function for each bin.

Parameters:
channelname The name of the channel.
processname The name of the process.
funccont A vector of pointers of TF1 functions.
nbins The number of bins used for the histogram.
efficiency The efficiency of this process in this channel.
Returns:
An error code.
See also:
SetTemplate(const char * channelname, const char * processname, TH1D hist, double efficiency = 1.)

Definition at line 162 of file BCMTF.cxx.

{
   // get channel index
   int channelindex = GetChannelIndex(channelname);

   // check if channel exists
   if (channelindex < 0) {
      BCLog::OutWarning("BCMultitemplateFitter::SetTemplate() : Channel does not exist.");
         return -1;
   }

   // get process index
   int processindex = GetProcessIndex(processname);

   // check if process exists
   if (processindex < 0) {
      BCLog::OutWarning("BCMultitemplateFitter::SetTemplate() : Process does not exist.");
         return -1;
   }

   // get channel
   BCMTFChannel * channel = GetChannel(channelindex);

   // get template
   BCMTFTemplate * bctemplate = channel->GetTemplate(processindex);

   // set histogram
   bctemplate->SetFunctionContainer(funccont, nbins);

   // set efficiency
   bctemplate->SetEfficiency(efficiency);

   // no error
   return 1;
}

Member Data Documentation

std::vector<BCMTFChannel *> BCMTF::fChannelContainer [private]

A container of channels.

Definition at line 382 of file BCMTF.h.

std::vector<TF1 *> BCMTF::fExpectationFunctionContainer [private]

A container of functions for the expectation.

Definition at line 420 of file BCMTF.h.

bool BCMTF::fFlagEfficiencyConstraint [private]

A flag for the efficiency constraint: efficiency within 0 and 1 (true) or not (false).

Definition at line 416 of file BCMTF.h.

int BCMTF::fNChannels [private]

The number of channels.

Definition at line 394 of file BCMTF.h.

int BCMTF::fNProcesses [private]

The number of processes.

Definition at line 398 of file BCMTF.h.

int BCMTF::fNSystematics [private]

The number of systematics uncertainties.

Definition at line 402 of file BCMTF.h.

std::vector<BCMTFProcess *> BCMTF::fProcessContainer [private]

A container of processes.

Definition at line 386 of file BCMTF.h.

std::vector<int> BCMTF::fProcessParIndexContainer [private]

A container of parameter indeces for the process normalization.

Definition at line 407 of file BCMTF.h.

std::vector<BCMTFSystematic *> BCMTF::fSystematicContainer [private]

A container of sources of systematic uncertainty.

Definition at line 390 of file BCMTF.h.

std::vector<int> BCMTF::fSystematicParIndexContainer [private]

A container of parameter indeces for the systematics.

Definition at line 411 of file BCMTF.h.

The documentation for this class was generated from the following files:
Generated by  doxygen 1.7.1