A class for doing ensemble tests. More...

#include <BCTemplateEnsembleTest.h>

Collaboration diagram for BCTemplateEnsembleTest:

Public Member Functions
	BCTemplateEnsembleTest ()
int	PerformEnsembleTest (TTree *tree=0)
int	PrepareTree ()
void	PrintPulls (const char *filename)
void	SetEnsembleExpectation (double expectation)
int	SetEnsembleTemplate (TH1D hist)
void	SetFlagMCMC (bool flag)
void	SetNEnsembles (int n)
void	SetTemplateFitter (BCTemplateFitter *model)
void	SetTemplateParameters (std::vector< double > parameters)
int	Write (const char *filename)
	~BCTemplateEnsembleTest ()
Protected Attributes
int	fEnsembleCounter
double	fEnsembleExpectation
TH1D	fEnsembleTemplate
TFile *	fFile
bool	fFlagMCMC
int	fNEnsembles
double	fOutChi2Global
double	fOutChi2Marg
double	fOutChi2ProbGlobal
double	fOutChi2ProbMarg
double	fOutKSProb
int	fOutNDF
int	fOutNEvents
std::vector< double >	fOutParErrorDownGlobal
std::vector< double >	fOutParErrorDownMarg
std::vector< double >	fOutParErrorUpGlobal
std::vector< double >	fOutParErrorUpMarg
std::vector< double >	fOutParMeanMarg
std::vector< double >	fOutParMedianMarg
std::vector< double >	fOutParModeGlobal
std::vector< double >	fOutParModeMarg
std::vector< double >	fOutParPullGlobal
std::vector< double >	fOutParPullMarg
std::vector< double >	fOutParQuantile10Marg
std::vector< double >	fOutParQuantile5Marg
std::vector< double >	fOutParQuantile90Marg
std::vector< double >	fOutParQuantile95Marg
std::vector< double >	fOutParRMSMarg
std::vector< double >	fOutParValue
double	fOutPValue
std::vector< double >	fOutRatioErrorDownMarg
std::vector< double >	fOutRatioErrorUpMarg
std::vector< double >	fOutRatioMeanMarg
std::vector< double >	fOutRatioMedianMarg
std::vector< double >	fOutRatioModeMarg
std::vector< double >	fOutRatioQuantile10Marg
std::vector< double >	fOutRatioQuantile5Marg
std::vector< double >	fOutRatioQuantile90Marg
std::vector< double >	fOutRatioQuantile95Marg
std::vector< double >	fOutRatioRMSMarg
TRandom3 *	fRandom
BCTemplateFitter *	fTemplateFitter
std::vector< double >	fTemplateParameters
TTree *	fTree
Private Member Functions
TH1D *	BuildEnsemble ()

Detailed Description

A class for doing ensemble tests.

This class can be used for ensemble tests using the StackTool. The fitting can be done with Minuit or with Markov Chains.

Author:: Andrea Knue; Daniel Kollar; Kevin Kroeninger

Date:: 10.04.2010

Definition at line 36 of file BCTemplateEnsembleTest.h.

Constructor & Destructor Documentation

BCTemplateEnsembleTest::BCTemplateEnsembleTest ( )

The constructor.

Definition at line 20 of file BCTemplateEnsembleTest.cxx.

   : fTemplateFitter(0)
   , fFile(0)
   , fTree(0)
   , fEnsembleCounter(0)
   , fEnsembleExpectation(0)
   , fNEnsembles(0)
   , fFlagMCMC(false)
{
   fRandom = new TRandom3(0);
}

BCTemplateEnsembleTest::~BCTemplateEnsembleTest ( )

The destructor.

Definition at line 33 of file BCTemplateEnsembleTest.cxx.

{
   if (fRandom)
      delete fRandom;
}

Member Function Documentation

TH1D * BCTemplateEnsembleTest::BuildEnsemble ( ) [private]

Create a new ensemble.

Returns:: A histogram with the new ensemble.

Definition at line 187 of file BCTemplateEnsembleTest.cxx.

{
   // get histogram parameters
  int nbins   = fTemplateFitter->GetData().GetNbinsX();

   // create new ensemble
   TH1D* ensemble = new TH1D(fTemplateFitter->GetData());

   // increase ensemble counter
   fEnsembleCounter++;

   // get new parameter if needed
   std::vector<double> parameters = fTemplateParameters;

   // loop over bins and fill them
  for(int ibin = 1; ibin <= nbins; ++ibin){
      double nexp = fTemplateFitter->Expectation(ibin, parameters);
      double nobs = gRandom->Poisson(nexp);

      // set the bin content
    ensemble->SetBinContent(ibin, nobs);
  }

   // return the ensemble histogram
  return ensemble;
}

int BCTemplateEnsembleTest::PerformEnsembleTest ( TTree * tree = 0 )

A function to perform an ensemble test for each data set in the container.

Parameters:

tree,:

A tree from which the parameters are taken, e.g., the poesterior

Definition at line 62 of file BCTemplateEnsembleTest.cxx.

{
   // set log level to nothing
   BCLog::LogLevel ll = BCLog::GetLogLevelScreen();
   BCLog::SetLogLevel(BCLog::nothing);

   // initialize template fitter
   fTemplateFitter->Initialize();

   // Prepare tree
   PrepareTree();

   // get number of parameters
   int npar = fTemplateFitter->GetNParameters();

   // connect tree
   if (tree) {
      fTemplateParameters = std::vector<double>(npar);
      for (int i = 0; i < npar; ++i) {
         tree->SetBranchAddress(Form("Parameter%i", i), &(fTemplateParameters[i]));
      }
   }

   // loop over ensembles
   for(int j = 0; j < fNEnsembles; j++){

      // print status
      if ((j+1) % 100 == 0 && j > 0)
         std::cout << "Fraction of ensembles analyzed: " << double(j+1) / double(fNEnsembles) * 100 << "%" << std::endl;

      // get parameters from tree
      if (tree) {
         int index = fRandom->Uniform(tree->GetEntries());
         tree->GetEntry(index);
      }

      // create new ensemble
      TH1D* ensemble = BuildEnsemble();

      // set ensemble as new data set
      fTemplateFitter->SetData(*ensemble);

      // find mode
      fTemplateFitter->FindMode();

      // perform MCMC
      if(fFlagMCMC) {
         fTemplateFitter->MCMCInitialize();

         fTemplateFitter->MarginalizeAll();

         // find mode with MCMC best fit
         fTemplateFitter->FindMode(fTemplateFitter->GetBestFitParameters());

         // loop over parameters and set tree variables
         for (int i = 0; i < npar; ++i) {
            BCH1D* hist = fTemplateFitter->GetMarginalized(fTemplateFitter->GetParameter(i));
            fOutParModeMarg[i]       = hist->GetMode();
            fOutParMedianMarg[i]     = hist->GetMedian();
            fOutParMeanMarg[i]       = hist->GetMean();
            fOutParRMSMarg[i]        = hist->GetRMS();
            fOutParErrorUpMarg[i]    = hist->GetQuantile(0.84)-hist->GetMode();
            fOutParErrorDownMarg[i]  = hist->GetMode()-hist->GetQuantile(0.16);
            fOutParQuantile5Marg[i]  = hist->GetQuantile(0.05);
            fOutParQuantile10Marg[i] = hist->GetQuantile(0.10);
            fOutParQuantile90Marg[i] = hist->GetQuantile(0.90);
            fOutParQuantile95Marg[i] = hist->GetQuantile(0.95);
            fOutParPullMarg[i]       = -1;
            double error = 0.5 * (fOutParErrorUpMarg.at(i) + fOutParErrorDownMarg.at(i));
            if (error > 0)
               fOutParPullMarg[i]     = (fOutParModeMarg.at(i) - fTemplateParameters.at(i)) / ( error );
         }
         fOutChi2Marg             = fTemplateFitter->CalculateChi2( fTemplateFitter->GetBestFitParametersMarginalized() );
         fOutChi2ProbMarg         = fTemplateFitter->CalculateChi2Prob(fTemplateFitter->GetBestFitParametersMarginalized());
      }

      if (fFlagMCMC) {
         int nratios = fTemplateFitter->GetNRatios();
         for (int i = 0; i < nratios; ++i) {
            TH1D histtemp = fTemplateFitter->GetHistRatio1D(i);
            BCH1D * hist = new BCH1D( &histtemp );
            fOutRatioModeMarg[i]       = hist->GetMode();
            fOutRatioMedianMarg[i]     = hist->GetMedian();
            fOutRatioMeanMarg[i]       = hist->GetMean();
            fOutRatioRMSMarg[i]        = hist->GetRMS();
            fOutRatioErrorUpMarg[i]    = hist->GetQuantile(0.84)-hist->GetMode();
            fOutRatioErrorDownMarg[i]  = hist->GetMode()-hist->GetQuantile(0.16);
            fOutRatioQuantile5Marg[i]  = hist->GetQuantile(0.05);
            fOutRatioQuantile10Marg[i] = hist->GetQuantile(0.10);
            fOutRatioQuantile90Marg[i] = hist->GetQuantile(0.90);
            fOutRatioQuantile95Marg[i] = hist->GetQuantile(0.95);
         }
      }

      // set tree variables
      fOutParValue           = fTemplateParameters;
      fOutParModeGlobal      = fTemplateFitter->GetBestFitParameters();
      fOutParErrorUpGlobal   = fTemplateFitter->GetBestFitParameterErrors();
      fOutParErrorDownGlobal = fTemplateFitter->GetBestFitParameterErrors();
      fOutChi2Global         = fTemplateFitter->CalculateChi2( fTemplateFitter->GetBestFitParameters() );
      fOutNDF                = fTemplateFitter->GetNDF();
      fOutChi2ProbGlobal     = fTemplateFitter->CalculateChi2Prob(fTemplateFitter->GetBestFitParameters());
      fOutKSProb             = fTemplateFitter->CalculateKSProb();
      fOutPValue             = fTemplateFitter->CalculatePValue();
      fOutNEvents            = int(fTemplateFitter->GetData().Integral());

      for (int i = 0; i < npar; ++i) {
         fOutParPullGlobal[i] = -1;
         double error = 0.5 * (fOutParErrorUpGlobal.at(i) + fOutParErrorDownGlobal.at(i));
         if (error > 0)
            fOutParPullGlobal[i] = (fOutParModeGlobal.at(i) - fTemplateParameters.at(i)) / ( error );
      }

      // fill the tree
      fTree->Fill();
   }

   // reset log level
   BCLog::SetLogLevel(ll);

   // no error
   return 1;
}

int BCTemplateEnsembleTest::PrepareTree ( )

Prepare tree.

Definition at line 234 of file BCTemplateEnsembleTest.cxx.

{
   // delete old tree if necessary
   if (fTree)
      delete fTree;

   // create new tree
   fTree = new TTree("fTree", "fTree");

   // get number of parameters and ratios
   int npar = fTemplateFitter->GetNParameters();
   int nratios = fTemplateFitter->GetNRatios();

   // initialize variables
   fOutParValue.assign(npar, 0);
   fOutParModeGlobal.assign(npar, 0);
   fOutParErrorUpGlobal.assign(npar, 0);
   fOutParErrorDownGlobal.assign(npar, 0);
   fOutParPullGlobal.assign(npar, 0);
   fOutParModeMarg.assign(npar, 0);
   fOutParMeanMarg.assign(npar, 0);
   fOutParMedianMarg.assign(npar, 0);
   fOutParRMSMarg.assign(npar, 0);
   fOutParErrorUpMarg.assign(npar, 0);
   fOutParErrorDownMarg.assign(npar, 0);
   fOutParPullMarg.assign(npar, 0);
   fOutParQuantile5Marg.assign(npar, 0);
   fOutParQuantile10Marg.assign(npar, 0);
   fOutParQuantile90Marg.assign(npar, 0);
   fOutParQuantile95Marg.assign(npar, 0);
   fOutRatioModeMarg.assign(nratios, 0);
   fOutRatioMeanMarg.assign(nratios, 0);
   fOutRatioMedianMarg.assign(nratios, 0);
   fOutRatioRMSMarg.assign(nratios, 0);
   fOutRatioErrorUpMarg.assign(nratios, 0);
   fOutRatioErrorDownMarg.assign(nratios, 0);
   fOutRatioQuantile5Marg.assign(nratios, 0);
   fOutRatioQuantile10Marg.assign(nratios, 0);
   fOutRatioQuantile90Marg.assign(nratios, 0);
   fOutRatioQuantile95Marg.assign(nratios, 0);

   fTree->Branch("chi2_global", &fOutChi2Global,     "chi2 (global)/D");
   fTree->Branch("ndf",      &fOutNDF,      "ndf/I");
   fTree->Branch("chi2prob_global", &fOutChi2ProbGlobal, "chi2 prob probability (global)/D");
   fTree->Branch("KSprob",   &fOutKSProb,   "KS probability/D");
   fTree->Branch("pvalue",   &fOutPValue,   "p-value/D");
   fTree->Branch("nevents",  &fOutNEvents,  "n events/I");

   for (int i = 0; i < npar; ++i) {
      // add branches
      fTree->Branch(Form("parameter_%i", i),                 &fOutParValue[i],      Form("parameter_%i/D", i));
      fTree->Branch(Form("par_global_mode_par_%i", i),       &fOutParModeGlobal[i],      Form("par_global_Mode_par_%i/D", i));
      fTree->Branch(Form("par_global_error_up_par_%i", i),   &fOutParErrorUpGlobal[i],   Form("par_global_error_up_par_%i/D", i));
      fTree->Branch(Form("par_global_error_down_par_%i", i), &fOutParErrorDownGlobal[i], Form("par_global_error_down_par_%i/D", i));
      fTree->Branch(Form("par_global_pull_par_%i", i),       &fOutParPullGlobal[i],      Form("par_global_pull_par_%i/D", i));

      if(fFlagMCMC) {
         fTree->Branch(Form("par_marg_mode_par_%i", i),       &fOutParModeMarg[i],       Form("par_marg_mode_par_%i/D", i));
         fTree->Branch(Form("par_marg_mean_par_%i", i),       &fOutParMeanMarg[i],       Form("par_marg_mean_par_%i/D", i));
         fTree->Branch(Form("par_marg_median_par_%i", i),     &fOutParMedianMarg[i],     Form("par_marg_median_par_%i/D", i));
         fTree->Branch(Form("par_marg_rms_par_%i", i),        &fOutParRMSMarg[i],        Form("par_marg_rms_par_%i/D", i));
         fTree->Branch(Form("par_marg_error_up_par_%i", i),   &fOutParErrorUpMarg[i],    Form("par_marg_ErrorUp_par_%i/D", i));
         fTree->Branch(Form("par_marg_error_down_par_%i", i), &fOutParErrorDownMarg[i],  Form("par_marg_error_down_par_%i/D", i));
         fTree->Branch(Form("par_marg_pull_par_%i", i),       &fOutParPullMarg[i],       Form("par_marg_pull_par_%i/D", i));
         fTree->Branch(Form("par_marg_quantile5_par_%i", i),  &fOutParQuantile5Marg[i],  Form("par_marg_Quantile5_par_%i/D", i));
         fTree->Branch(Form("par_marg_quantile10_par_%i", i), &fOutParQuantile10Marg[i], Form("par_marg_Quantile10_par_%i/D", i));
         fTree->Branch(Form("par_marg_quantile90_par_%i", i), &fOutParQuantile90Marg[i], Form("par_marg_Quantile90_par_%i/D", i));
         fTree->Branch(Form("par_marg_quantile95_par_%i", i), &fOutParQuantile95Marg[i], Form("par_marg_Quantile95_par_%i/D", i));
      }
   }

   if (fFlagMCMC) {
      fTree->Branch("chi2_marg", &fOutChi2Marg,     "chi2 (marginalized)/D");
      fTree->Branch("chi2prob_marg", &fOutChi2ProbMarg, "chi2 prob probability (marginalized)/D");

      for (int i = 0; i < nratios; ++i) {
         fTree->Branch(Form("ratio_marg_mode_ratio_%i", i),       &fOutRatioModeMarg[i],       Form("ratio_marg_mode_ratio_%i/D", i));
         fTree->Branch(Form("ratio_marg_mean_ratio_%i", i),       &fOutRatioMeanMarg[i],       Form("ratio_marg_mean_ratio_%i/D", i));
         fTree->Branch(Form("ratio_marg_median_ratio_%i", i),     &fOutRatioMedianMarg[i],     Form("ratio_marg_median_ratio_%i/D", i));
         fTree->Branch(Form("ratio_marg_rms_ratio_%i", i),        &fOutRatioRMSMarg[i],        Form("ratio_marg_rms_ratio_%i/D", i));
         fTree->Branch(Form("ratio_marg_error_up_ratio_%i", i),   &fOutRatioErrorUpMarg[i],    Form("ratio_marg_ErrorUp_ratio_%i/D", i));
         fTree->Branch(Form("ratio_marg_error_down_ratio_%i", i), &fOutRatioErrorDownMarg[i],  Form("ratio_marg_error_down_ratio_%i/D", i));
         fTree->Branch(Form("ratio_marg_quantile5_ratio_%i", i),  &fOutRatioQuantile5Marg[i],  Form("ratio_marg_Quantile5_ratio_%i/D", i));
         fTree->Branch(Form("ratio_marg_quantile10_ratio_%i", i), &fOutRatioQuantile10Marg[i], Form("ratio_marg_Quantile10_ratio_%i/D", i));
         fTree->Branch(Form("ratio_marg_quantile90_ratio_%i", i), &fOutRatioQuantile90Marg[i], Form("ratio_marg_Quantile90_ratio_%i/D", i));
         fTree->Branch(Form("ratio_marg_quantile95_ratio_%i", i), &fOutRatioQuantile95Marg[i], Form("ratio_marg_Quantile95_ratio_%i/D", i));
      }
   }

   // no error
   return 1;
}

void BCTemplateEnsembleTest::PrintPulls ( const char * filename )

Print pulls

Definition at line 329 of file BCTemplateEnsembleTest.cxx.

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

   // create canvas and prepare postscript
   TCanvas * canvas = new TCanvas();

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

   // get number of parameters
   int npar = fTemplateFitter->GetNParameters();

   // create histogram container
   std::vector<TH1D*> histcontainer = std::vector<TH1D*>(0);

   // loop over all parameters
   for (int j = 0; j < npar; ++j) { 
      TH1D* hist = new TH1D("", "Pull;N", 11, -5.5, 5.5);
      histcontainer.push_back(hist);
   }

   // loop over all ensembles
   for (int i = 0; i < fNEnsembles; ++i) {

      // get ensemble
      fTree->GetEntry(i);

      // loop over all parameters
      for (int j = 0; j < npar; ++j) { 
         histcontainer[j]->Fill(fOutParPullGlobal[j]);
      }
   }

   // loop over all parameters
   for (int j = 0; j < npar; ++j) { 
      // update post script
      canvas->Update();
      ps->NewPage();
      canvas->cd();

      canvas->cd();
      histcontainer.at(j)->Draw();
   }
   // close ps
   canvas->Update();
   ps->Close();

   // free memory
   delete ps;
   delete canvas;

}

void BCTemplateEnsembleTest::SetEnsembleExpectation ( double expectation ) [inline]

A function to define the number of events per ensemble.

Definition at line 79 of file BCTemplateEnsembleTest.h.

      { fEnsembleExpectation = expectation; };

int BCTemplateEnsembleTest::SetEnsembleTemplate ( TH1D hist )

Set the template used to generate the ensembles.

Definition at line 40 of file BCTemplateEnsembleTest.cxx.

{
   // calculate integral
   double integral = hist.Integral();

   // check if integral is ok
   if (integral <= 0) {
      std::cout << "Template not valid. Integral is lower or equal to 0." << std::endl;
      return 0;
   }

   // set template
   fEnsembleTemplate = hist;

   // scale template
   fEnsembleTemplate.Scale(1.0/integral);

   // no error
   return 1;
};

void BCTemplateEnsembleTest::SetFlagMCMC ( bool flag ) [inline]

A function to set the MCMC flag.

Definition at line 91 of file BCTemplateEnsembleTest.h.

     { fFlagMCMC = flag; }

void BCTemplateEnsembleTest::SetNEnsembles ( int n ) [inline]

A function to define the number of ensembles per data set.

Definition at line 85 of file BCTemplateEnsembleTest.h.

      { fNEnsembles = n; };

void BCTemplateEnsembleTest::SetTemplateFitter ( BCTemplateFitter * model ) [inline]

Set the BCTemplateFitter used to analyze the ensembles.

Definition at line 58 of file BCTemplateEnsembleTest.h.

      { fTemplateFitter = model; };

void BCTemplateEnsembleTest::SetTemplateParameters ( std::vector< double > parameters ) [inline]

Set the parameters used for generating ensembles using the current template fitter model.

Definition at line 65 of file BCTemplateEnsembleTest.h.

         { fTemplateParameters = parameters; };

int BCTemplateEnsembleTest::Write ( const char * filename )

Write tree to file

Definition at line 215 of file BCTemplateEnsembleTest.cxx.

{
   // open file
   fFile = new TFile(filename, "RECREATE");

   // write tree
   fTree->Write();

   // close file
   fFile->Close();

   // free memory
   delete fFile;

   // no error
   return 1;
}