BCMTFAnalysisFacility.cxx

/*
 * Copyright (C) 2007-2018, the BAT core developer team
 * All rights reserved.
 *
 * For the licensing terms see doc/COPYING.
 * For documentation see http://mpp.mpg.de/bat
 */

// ---------------------------------------------------------

#include "BCMTFAnalysisFacility.h"

#include "BCMTF.h"
#include "BCMTFChannel.h"
#include "BCMTFComparisonTool.h"
#include "BCMTFSystematic.h"
#include "BCMTFTemplate.h"

#include "../../BAT/BCLog.h"
#include "../../BAT/BCH1D.h"
#include "../../BAT/BCParameter.h"

#include <TCanvas.h>
#include <TFile.h>
#include <TH1D.h>
#include <TRandom3.h>
#include <TROOT.h>
#include <TTree.h>

#include <sys/stat.h>
#include <unistd.h>
#include <iostream>
#include <stdexcept>

namespace
{
void cd(const std::string& dir)
{
    // shell conventions: 0=success
    if (chdir(dir.data()))
        throw std::runtime_error(std::string("Cannot change directory to ") + dir);
}
}

// ---------------------------------------------------------
BCMTFAnalysisFacility::BCMTFAnalysisFacility(BCMTF* mtf)
    : fRandom(new TRandom3(0))
    , fFlagMarginalize(false)
    , fLogLevel(BCLog::nothing)
{
    fMTF = mtf;
    BCLog::OutDetail("Prepared Analysis Facility for MTF model \'" + mtf->GetName() + "\'");
}

// ---------------------------------------------------------
BCMTFAnalysisFacility::~BCMTFAnalysisFacility()
{
}

// ---------------------------------------------------------
std::vector<TH1D> BCMTFAnalysisFacility::BuildEnsemble(const std::vector<double>& parameters, const std::string& options)
{
    // option flags
    bool flag_data = false;

    // check content of options string
    if (options.find("data") < options.size()) {
        flag_data = true;
    }

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

    // create vector of histograms
    std::vector<TH1D> histograms;

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

        // get channel
        BCMTFChannel* channel = fMTF->GetChannel(ichannel);

        // create new histogram
        TH1D hist( *(channel->GetData()->GetHistogram()) );

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

        // loop over all bins
        for (int ibin = 1; ibin <= nbins; ++ibin) {
            if (!flag_data) {
                double expectation = fMTF->Expectation(ichannel, ibin, parameters);
                double observation = fRandom->Poisson(expectation);

                hist.SetBinContent(ibin, observation);
            }
        }

        // add histogram
        histograms.push_back(hist);
    }

    // return histograms
    return histograms;
}

// ---------------------------------------------------------
TTree* BCMTFAnalysisFacility::BuildEnsembles(TTree* tree, int nensembles, const std::string& options)
{
    // get number of channels
    int nchannels = fMTF->GetNChannels();

    BCLog::OutDetail(Form("MTF Building %d ensembles for %d channels.", nensembles, nchannels));

    // get number of parameters
    int nparameters = fMTF->GetNParameters();

    // create tree variables for input tree
    std::vector<double> parameters(nparameters);

    // set branch addresses
    for (int i = 0; i < nparameters; ++i) {
        tree->SetBranchAddress(Form("Parameter%i", i), &(parameters[i]));
    }

    // create tree
    TTree* tree_out = new TTree("ensembles", "ensembles");

    // create matrix of number of bins
    std::vector< std::vector<double> > nbins_matrix;

    // prepare the tree variables
    // loop over channels
    for (int ichannel = 0; ichannel < nchannels; ++ichannel) {
        // get channel
        BCMTFChannel* channel = fMTF->GetChannel(ichannel);

        // get number of bins
        int nbins = channel->GetData()->GetHistogram()->GetNbinsX();

        // create new matrix row
        std::vector<double> nbins_column(nbins);

        // add matrix
        nbins_matrix.push_back(nbins_column);
    }

    std::vector<double> in_parameters(nparameters);

    // create branches
    // loop over channels
    for (int ichannel = 0; ichannel < nchannels; ++ichannel) {
        // get channel
        BCMTFChannel* channel = fMTF->GetChannel(ichannel);

        // get number of bins
        int nbins = channel->GetData()->GetHistogram()->GetNbinsX();

        // loop over bins
        for (int ibin = 1; ibin <= nbins; ++ibin) {
            // create branches
            tree_out->Branch(Form("channel_%i_bin_%i", ichannel, ibin),
                             &(nbins_matrix[ichannel])[ibin - 1], "n/D");
        }
    }

    for (int i = 0; i < nparameters; ++i) {
        tree_out->Branch(Form("parameter_%i", i), &in_parameters[i], Form("parameter_%i/D", i));
    }

    // create vector of histograms
    std::vector<TH1D> histograms;

    // loop over ensembles
    for (int iensemble = 0; iensemble < nensembles; ++iensemble) {
        // get random event from tree
        int index = (int) fRandom->Uniform(tree->GetEntries());
        tree->GetEntry(index);

        // create ensembles
        histograms = BuildEnsemble(parameters, options);

        // copy information from histograms into tree variables
        // loop over channels
        for (int ichannel = 0; ichannel < nchannels; ++ichannel) {
            // get channel
            BCMTFChannel* channel = fMTF->GetChannel(ichannel);

            // get number of bins
            int nbins = channel->GetData()->GetHistogram()->GetNbinsX();

            // loop over bins
            for (int ibin = 1; ibin <= nbins; ++ibin) {
                // fill tree variable
                (nbins_matrix[ichannel])[ibin - 1] = histograms.at(ichannel).GetBinContent(ibin);
            }
        }

        // copy parameter information
        for (int i = 0; i < nparameters; ++i) {
            in_parameters[i] = parameters.at(i);
        }

        // fill tree
        tree_out->Fill();
    }

    // return tree
    return tree_out;
}

// ---------------------------------------------------------
TTree* BCMTFAnalysisFacility::BuildEnsembles(const std::vector<double>& parameters, int nensembles, const std::string& options)
{
    // get number of channels
    int nchannels = fMTF->GetNChannels();

    BCLog::OutDetail(Form("MTF Building %d ensambles for %d channels.", nensembles, nchannels));

    // create tree
    TTree* tree = new TTree("ensembles", "ensembles");

    // create matrix of number of bins
    std::vector< std::vector<double> > nbins_matrix;

    // prepare the tree variables
    // loop over channels
    for (int ichannel = 0; ichannel < nchannels; ++ichannel) {
        // get channel
        BCMTFChannel* channel = fMTF->GetChannel(ichannel);

        // get number of bins
        int nbins = channel->GetData()->GetHistogram()->GetNbinsX();

        // create new matrix row
        std::vector<double> nbins_column(nbins);

        // add matrix
        nbins_matrix.push_back(nbins_column);
    }

    // get number of parameters
    int nparameters = fMTF->GetNParameters();

    std::vector<double> in_parameters(nparameters);

    // create branches
    // loop over channels
    for (int ichannel = 0; ichannel < nchannels; ++ichannel) {
        // get channel
        BCMTFChannel* channel = fMTF->GetChannel(ichannel);

        // get number of bins
        int nbins = channel->GetData()->GetHistogram()->GetNbinsX();

        // loop over bins
        for (int ibin = 1; ibin <= nbins; ++ibin) {
            // create branches
            tree->Branch(Form("channel_%i_bin_%i", ichannel, ibin),
                         &(nbins_matrix[ichannel])[ibin - 1], "n/D");
        }
    }

    for (int i = 0; i < nparameters; ++i) {
        tree->Branch(Form("parameter_%i", i), &in_parameters[i], Form("parameter_%i/D", i));
    }

    // create vector of histograms
    std::vector<TH1D> histograms;

    // loop over ensembles
    for (int iensemble = 0; iensemble < nensembles; ++iensemble) {
        // create ensembles
        histograms = BuildEnsemble(parameters, options);

        // copy information from histograms into tree variables
        // loop over channels
        for (int ichannel = 0; ichannel < nchannels; ++ichannel) {
            // get channel
            BCMTFChannel* channel = fMTF->GetChannel(ichannel);

            // get number of bins
            int nbins = channel->GetData()->GetHistogram()->GetNbinsX();

            // loop over bins
            for (int ibin = 1; ibin <= nbins; ++ibin) {
                // fill tree variable
                (nbins_matrix[ichannel])[ibin - 1] = histograms.at(ichannel).GetBinContent(ibin);
            }
        }

        // copy parameter information
        for (int i = 0; i < nparameters; ++i) {
            if (parameters.size() > 0)
                in_parameters[i] = parameters.at(i);
            else
                in_parameters[i] = 0;
        }

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

    // return tree
    return tree;
}

// ---------------------------------------------------------
TTree* BCMTFAnalysisFacility::PerformEnsembleTest(const std::vector<double>& parameters, int nensembles, const std::string& options)
{
    // create new tree
    TTree* tree = 0;

    // create ensembles
    tree = BuildEnsembles(parameters, nensembles, options);

    // perform ensemble test
    return PerformEnsembleTest(tree, nensembles, 0, options);
}

// ---------------------------------------------------------
TTree* BCMTFAnalysisFacility::PerformEnsembleTest(TTree* tree, int nensembles, int start, const std::string& options)
{
    BCLog::OutSummary("Running ensemble test.");
    if (fFlagMarginalize) {
        BCLog::OutSummary("Fit for each ensemble is going to be run using MCMC. It can take a while.");
    }

    // option flags
    bool flag_mc = false;

    // check content of options string
    if (options.find("MC") < options.size()) {
        flag_mc = true;
    }

    // set log level
    // It would be better to set the level for the screen and for the file
    // separately. Perhaps in the future.
    BCLog::LogLevel lls = BCLog::GetLogLevelScreen();
    BCLog::LogLevel llf = BCLog::GetLogLevelFile();
    if (fLogLevel == BCLog::nothing) {
        BCLog::OutSummary("No log messages for the ensemble fits are going to be printed.");
        BCLog::SetLogLevel(fLogLevel);
    } else if (fLogLevel != lls) {
        BCLog::OutSummary("The log level for the ensemble test is set to '" + BCLog::ToString(fLogLevel) + "'.");
        BCLog::SetLogLevel(fLogLevel);
    }

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

    // define set of histograms for the original data sets
    std::vector<TH1D*> histograms_data(nchannels);

    // create matrix of number of bins
    std::vector< std::vector<double> > nbins_matrix;

    // prepare the tree
    // loop over channels
    for (int ichannel = 0; ichannel < nchannels; ++ichannel) {
        // get channel
        BCMTFChannel* channel = fMTF->GetChannel(ichannel);

        // get number of bins
        int nbins = channel->GetData()->GetHistogram()->GetNbinsX();

        // create new matrix row
        std::vector<double> nbins_column(nbins);

        // add matrix
        nbins_matrix.push_back(nbins_column);
    }

    // loop over channels
    for (int ichannel = 0; ichannel < nchannels; ++ichannel) {
        // get channel
        BCMTFChannel* channel = fMTF->GetChannel(ichannel);

        // get number of bins
        int nbins = channel->GetData()->GetHistogram()->GetNbinsX();

        // loop over bins
        for (int ibin = 1; ibin <= nbins; ++ibin) {
            // create branches
            tree->SetBranchAddress(Form("channel_%i_bin_%i", ichannel, ibin),
                                   &(nbins_matrix[ichannel])[ibin - 1]);
        }
    }

    // get number of parameters
    int nparameters = fMTF->GetNParameters();

    // define tree variables
    std::vector<double> out_parameters(nparameters);

    for (int i = 0; i < nparameters; ++i) {
        tree->SetBranchAddress(Form("parameter_%i", i), &out_parameters[i]);
    }

    // copy the original data sets
    // loop over channels
    for (int ichannel = 0; ichannel < nchannels; ++ichannel) {
        // get channel
        BCMTFChannel* channel = fMTF->GetChannel(ichannel);

        // set data pointer
        histograms_data[ichannel] = channel->GetData()->GetHistogram();
    }

    // create output tree
    TTree* tree_out = new TTree("ensemble_test", "ensemble test");

    // define tree variables
    std::vector<double> out_mode_global(nparameters);
    std::vector<double> out_std_global(nparameters);
    std::vector<double> out_mode_marginalized(nparameters);
    std::vector<double> out_mean_marginalized(nparameters);
    std::vector<double> out_median_marginalized(nparameters);
    std::vector<double> out_5quantile_marginalized(nparameters);
    std::vector<double> out_10quantile_marginalized(nparameters);
    std::vector<double> out_16quantile_marginalized(nparameters);
    std::vector<double> out_84quantile_marginalized(nparameters);
    std::vector<double> out_90quantile_marginalized(nparameters);
    std::vector<double> out_95quantile_marginalized(nparameters);
    std::vector<double> out_std_marginalized(nparameters);
    std::vector<double> out_chi2_generated(nchannels);
    std::vector<double> out_chi2_mode(nchannels);
    std::vector<double> out_cash_generated(nchannels);
    std::vector<double> out_cash_mode(nchannels);
    std::vector<int> out_nevents(nchannels);
    double out_chi2_generated_total;
    double out_chi2_mode_total;
    double out_cash_generated_total;
    double out_cash_mode_total;
    int out_nevents_total;

    // create branches
    for (int i = 0; i < nparameters; ++i) {
        tree_out->Branch(Form("parameter_%i", i), &out_parameters[i], Form("parameter %i/D", i));
        tree_out->Branch(Form("mode_global_%i", i), &out_mode_global[i], Form("global mode of par. %i/D", i));
        tree_out->Branch(Form("std_global_%i", i), &out_std_global[i], Form("global std of par. %i/D", i));
        if (fFlagMarginalize) {
            tree_out->Branch(Form("mode_marginalized_%i", i), &out_mode_marginalized[i], Form("marginalized mode of par. %i/D", i));
            tree_out->Branch(Form("mean_marginalized_%i", i), &out_mean_marginalized[i], Form("marginalized mean of par. %i/D", i));
            tree_out->Branch(Form("median_marginalized_%i", i), &out_median_marginalized[i], Form("median of par. %i/D", i));
            tree_out->Branch(Form("5quantile_marginalized_%i", i), &out_5quantile_marginalized[i], Form("marginalized 5 per cent quantile of par. %i/D", i));
            tree_out->Branch(Form("10quantile_marginalized_%i", i), &out_10quantile_marginalized[i], Form("marginalized 10 per cent quantile of par. %i/D", i));
            tree_out->Branch(Form("16quantile_marginalized_%i", i), &out_16quantile_marginalized[i], Form("marginalized 16 per cent quantile of par. %i/D", i));
            tree_out->Branch(Form("84quantile_marginalized_%i", i), &out_84quantile_marginalized[i], Form("marginalized 84 per cent quantile of par. %i/D", i));
            tree_out->Branch(Form("90quantile_marginalized_%i", i), &out_90quantile_marginalized[i], Form("marginalized 90 per cent quantile of par. %i/D", i));
            tree_out->Branch(Form("95quantile_marginalized_%i", i), &out_95quantile_marginalized[i], Form("marginalized 95 per cent quantile of par. %i/D", i));
            tree_out->Branch(Form("std_marginalized_%i", i), &out_std_marginalized[i], Form("marginalized std of par. %i/D", i));
        }
    }
    for (int i = 0; i < nchannels; ++i) {
        tree_out->Branch(Form("chi2_generated_%i", i), &out_chi2_generated[i], Form("chi2 (generated par.) in channel %i/D", i));
        tree_out->Branch(Form("chi2_mode_%i", i), &out_chi2_mode[i], Form("chi2 (mode of par.) in channel %i/D", i));
        tree_out->Branch(Form("cash_generated_%i", i), &out_cash_generated[i], Form("cash statistic (generated par.) in channel %i/D", i));
        tree_out->Branch(Form("cash_mode_%i", i), &out_cash_mode[i], Form("cash statistic (mode of par.) in channel %i/D", i));
        tree_out->Branch(Form("nevents_%i", i), &out_nevents[i], Form("nevents in channel %i/I", i));
    }
    tree_out->Branch("chi2_generated_total", &out_chi2_generated_total, "chi2 (generated par.) in all channels/D");
    tree_out->Branch("chi2_mode_total", &out_chi2_mode_total, "chi2 (mode of par.) in all channels/D");
    tree_out->Branch("cash_generated_total", &out_cash_generated_total, "cash statistic (generated par.) in all channels/D");
    tree_out->Branch("cash_mode_total", &out_cash_mode_total, "cash statistic (mode of par.) in all channels/D");
    tree_out->Branch("nevents_total", &out_nevents_total, "total number of events/I");

    // define temporary vector of histogram for fluctated templates
    std::vector<TH1D*> histlist(0);

    // loop over ensembles
    for (int iensemble = 0; iensemble < nensembles; ++iensemble) {
        // print status
        if ((iensemble + 1) % 100 == 0 && iensemble > 0) {
            BCLog::SetLogLevel(lls, llf);
            int frac = int (double(iensemble + 1) / double(nensembles) * 100.);
            BCLog::OutDetail(Form("Fraction of ensembles analyzed: %i%%", frac));
            BCLog::SetLogLevel(fLogLevel);
        }

        // get next (commented out: random) event from tree
        //      int index = (int) fRandom->Uniform(tree->GetEntries());
        int index = iensemble + start;
        tree->GetEntry(index);

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

        // transform matrix into histograms
        histograms = MatrixToHistograms(nbins_matrix);

        // loop over channels and set data
        for (int ichannel = 0; ichannel < nchannels; ++ichannel) {
            // get channel
            BCMTFChannel* channel = fMTF->GetChannel(ichannel);

            // overwrite data
            if (iensemble == 0)
                channel->GetData()->SetHistogram(0);

            // set data histogram
            fMTF->SetData(channel->GetName(), histograms.at(ichannel));
        }

        // fluctuate templates if option "MC" is chosen
        if (flag_mc) {
            // get number of templates
            unsigned int ntemplates = fMTF->GetNProcesses();

            // loop over channels
            for (int ichannel = 0; ichannel < nchannels; ++ichannel) {
                // get channel
                BCMTFChannel* channel = fMTF->GetChannel(ichannel);

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

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

                    // replace by fluctuated histogram
                    if (temphist) {
                        TH1D* temphistfluc = new TH1D(channel->GetTemplate(i)->FluctuateHistogram(options, channel->GetTemplate(i)->GetOriginalNorm()));
                        channel->GetTemplate(i)->SetHistogram(temphistfluc, channel->GetTemplate(i)->GetNorm());
                    }
                }
            }
        }

        // work-around: force initialization
        fMTF->ResetResults();

        // check if marginalization should be run and perform analysis
        if (fFlagMarginalize) {
            BCLog::SetLogLevel(lls, llf);
            BCLog::OutDetail(Form("Running MCMC for ensemble %i", iensemble));
            BCLog::SetLogLevel(fLogLevel);

            // run mcmc
            fMTF->MarginalizeAll();

            // find mode
            fMTF->FindMode(fMTF->GetBestFitParameters());
        } else {
            // find mode
            fMTF->FindMode();
        }

        // fill tree variables
        out_mode_global = fMTF->GetBestFitParameters();
        out_std_global = fMTF->GetBestFitParameterErrors();
        out_chi2_generated_total = 0;
        out_chi2_mode_total = 0;
        out_cash_generated_total = 0;
        out_cash_mode_total = 0;
        out_nevents_total = 0;

        double quantile_probsums[7] = {5e-2, 10e-2, 16e-2, 50e-2, 84e-2, 90e-2, 95e-2};
        double quantile_values[7];

        for (int i = 0; i < nparameters; ++i) {
            if (fFlagMarginalize) {
                BCH1D hist = fMTF->GetMarginalized(i);
                if (!hist.Valid())
                    continue;

                out_mode_marginalized[i]       = hist.GetLocalMode(0);
                out_mean_marginalized[i]       = hist.GetHistogram()->GetMean();
                out_std_marginalized[i]        = hist.GetHistogram()->GetRMS();

                hist.GetHistogram()->GetQuantiles(6, quantile_values, quantile_probsums);
                out_5quantile_marginalized[i]  = quantile_values[0];
                out_10quantile_marginalized[i] = quantile_values[1];
                out_16quantile_marginalized[i] = quantile_values[2];
                out_median_marginalized[i]     = quantile_values[3];
                out_84quantile_marginalized[i] = quantile_values[4];
                out_90quantile_marginalized[i] = quantile_values[5];
                out_95quantile_marginalized[i] = quantile_values[6];
            }
        }

        for (int ichannel = 0; ichannel < nchannels; ++ichannel) {
            // get channel
            BCMTFChannel* channel = fMTF->GetChannel(ichannel);

            // get number of events
            out_nevents[ichannel] = (int) channel->GetData()->GetHistogram()->Integral();

            // calculate chi2
            out_chi2_generated[ichannel] = fMTF->CalculateChi2( ichannel, out_parameters );
            out_chi2_mode[ichannel] = fMTF->CalculateChi2( ichannel, fMTF->GetBestFitParameters() );

            // calculate cash statistic
            out_cash_generated[ichannel] = fMTF->CalculateCash( ichannel, out_parameters );
            out_cash_mode[ichannel] = fMTF->CalculateCash( ichannel, fMTF->GetBestFitParameters() );

            // increase the total number of events
            out_nevents_total += out_nevents[ichannel];

            // increase the total chi2
            out_chi2_generated_total += out_chi2_generated[ichannel];
            out_chi2_mode_total += out_chi2_mode[ichannel];

            // increase the total cash statistic
            out_cash_generated_total += out_cash_generated[ichannel];
            out_cash_mode_total += out_cash_mode[ichannel];
        }

        // fill tree
        tree_out->Fill();

        // but original template back if options "MC" is chosen
        if (flag_mc) {
            // get number of templates
            unsigned int ntemplates = fMTF->GetNProcesses();

            // loop over channels
            for (int ichannel = 0; ichannel < nchannels; ++ichannel) {
                // get channel
                BCMTFChannel* channel = fMTF->GetChannel(ichannel);

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

                    // get histogram
                    TH1D* temphist = histlist.at(ichannel * ntemplates + i);
                    temphist->Scale(channel->GetTemplate(i)->GetOriginalNorm() / temphist->Integral());

                    // replace by fluctuated histogram
                    TH1D* temphistfluc = channel->GetTemplate(i)->GetHistogram();
                    delete temphistfluc;
                    channel->GetTemplate(i)->SetHistogram(temphist, channel->GetTemplate(i)->GetNorm());
                }
            }
            histlist.clear();
        }
    }

    // put the original data back in place
    // loop over channels
    for (int ichannel = 0; ichannel < nchannels; ++ichannel) {
        // get channel
        BCMTFChannel* channel = fMTF->GetChannel(ichannel);

        // set data pointer
        channel->GetData()->SetHistogram(histograms_data.at(ichannel));
    }

    // reset log level
    BCLog::SetLogLevel(lls, llf);

    // work-around: force initialization
    fMTF->ResetResults();

    BCLog::OutSummary("Ensemble test ran successfully.");

    // return output tree
    return tree_out;
}

// ---------------------------------------------------------
std::vector<TH1D> BCMTFAnalysisFacility::MatrixToHistograms(const std::vector< std::vector<double> >& matrix)
{
    // create vector of histograms
    std::vector<TH1D> histograms;

    // get number of channels
    int nchannels = matrix.size();;

    // loop over channels
    for (int ichannel = 0; ichannel < nchannels; ++ichannel) {
        // get channel
        BCMTFChannel* channel = fMTF->GetChannel(ichannel);

        // get column
        std::vector<double> nbins_column = matrix[ichannel];

        // create new histogram
        TH1D hist( *(channel->GetData()->GetHistogram()) );

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

        // fill bin content
        for (int ibin = 1; ibin <= nbins; ++ibin) {
            hist.SetBinContent(ibin, nbins_column.at(ibin - 1));
        }

        // add histogram to container
        histograms.push_back(hist);
    }

    // return histograms
    return histograms;

}

// ---------------------------------------------------------
void BCMTFAnalysisFacility::PerformSingleChannelAnalyses(const std::string& dirname, const std::string& options)
{
    BCLog::OutSummary("Running single channel analysis in directory \'" + dirname + "\'");

    // ---- create new directory ---- //
    mkdir(dirname.data(), 0777);
    cd(dirname);

    // ---- check options ---- //

    bool flag_syst = true;
    bool flag_mcmc = true;

    if (options.find("nosyst") < options.size())
        flag_syst = false;

    if (options.find("mcmc") < options.size())
        flag_mcmc = true;

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

    // get number of systematics
    int nsystematics = fMTF->GetNSystematics();

    // container of flags for channels
    std::vector<bool> flag_channel(nchannels);

    // container of flags for systematics
    std::vector<bool> flag_systematic(nsystematics);

    // create new container of comparison tools
    std::vector<BCMTFComparisonTool*> ctc;
    std::vector<BCMTFComparisonTool*> ctc_mcmc;

    // get number of parameters
    int nparameters = fMTF->GetNParameters();

    // ---- add one comparison tool for each parameter ---- //

    // loop over all parameters
    for (int i = 0; i < nparameters; ++ i) {
        // create new comparison tool
        BCMTFComparisonTool* ct = new BCMTFComparisonTool(fMTF->GetParameter(i).GetName());
        BCMTFComparisonTool* ct_mcmc = new BCMTFComparisonTool(fMTF->GetParameter(i).GetName());

        // add comparison tool to container
        ctc.push_back(ct);
        ctc_mcmc.push_back(ct_mcmc);
    }

    // ---- switch on/off all systematics ---- //

    // loop over all systematics
    for (int isystematic = 0; isystematic < nsystematics; ++isystematic) {
        // get systematic
        BCMTFSystematic* systematic = fMTF->GetSystematic(isystematic);

        // remember old setting
        flag_systematic[isystematic] = systematic->GetFlagSystematicActive();

        // switch off systematic
        if (flag_systematic[isystematic])
            systematic->SetFlagSystematicActive(flag_syst);
    }

    // ---- switch on all channels ---- //

    // loop over all channels
    for (int ichannel = 0; ichannel < nchannels; ++ichannel) {
        // get channel
        BCMTFChannel* channel = fMTF->GetChannel(ichannel);

        // remember old setting
        flag_channel[ichannel] = channel->GetFlagChannelActive();

        // switch channel on/off
        channel->SetFlagChannelActive(true);
    }

    // ---- perform analysis for combination ---- //
    if (flag_mcmc) {
        // work-around: force initialization
        fMTF->ResetResults();

        // run mcmc
        fMTF->MarginalizeAll();

        // find mode
        fMTF->FindMode(fMTF->GetBestFitParameters());
    } else {
        // find mode
        fMTF->FindMode();
    }

    // print results
    if (flag_mcmc)
        fMTF->PrintAllMarginalized("marginalized_combined.pdf");
    fMTF->PrintFitSummary();

    // loop over all parameters
    for (int i = 0; i < nparameters; ++ i) {
        // get comparison tool
        BCMTFComparisonTool* ct = ctc.at(i);
        BCMTFComparisonTool* ct_mcmc = ctc_mcmc.at(i);

        ct->AddContribution("all channels", fMTF->GetBestFitParameters().at(i), fMTF->GetBestFitParameterErrors().at(i));
        if (flag_mcmc) {
            TH1* hist = fMTF->GetMarginalizedHistogram(i);
            if (hist)
                ct_mcmc->AddContribution("all channels", hist->GetMean(), hist->GetRMS());
        }
    }

    // ---- switch off all channels ----//

    // loop over all channels
    for (int ichannel = 0; ichannel < nchannels; ++ichannel) {
        // get channel
        BCMTFChannel* channel = fMTF->GetChannel(ichannel);

        // switch off channel
        channel->SetFlagChannelActive(false);
    }

    // ---- perform analysis on all channels separately ---- //

    // loop over all channels
    for (int ichannel = 0; ichannel < nchannels; ++ichannel) {
        // get channel
        BCMTFChannel* channel = fMTF->GetChannel(ichannel);

        // switch on one channel
        channel->SetFlagChannelActive(true);

        // perform analysis

        if (flag_mcmc) {
            // work-around: force initialization
            fMTF->ResetResults();

            // run mcmc
            fMTF->MarginalizeAll();

            // find mode
            fMTF->FindMode(fMTF->GetBestFitParameters());
        } else {
            // find mode
            fMTF->FindMode();
        }

        // print results
        if (flag_mcmc)
            fMTF->PrintAllMarginalized(Form("marginalized_channel_%i.pdf", ichannel));
        BCLog::OutSummary(Form("Summary for fitting channel %i", ichannel));
        fMTF->PrintFitSummary();

        // ---- update comparison tools ---- //

        // loop over all parameters
        for (int i = 0; i < nparameters; ++ i) {
            // get comparison tool
            BCMTFComparisonTool* ct = ctc.at(i);
            BCMTFComparisonTool* ct_mcmc = ctc_mcmc.at(i);

            ct->AddContribution(channel->GetName(), fMTF->GetBestFitParameters().at(i), fMTF->GetBestFitParameterErrors().at(i));
            if (flag_mcmc) {
                TH1* hist = fMTF->GetMarginalizedHistogram(i);
                if (hist)
                    ct_mcmc->AddContribution(channel->GetName(), hist->GetMean(), hist->GetRMS());
            }

            // switch off channel
            channel->SetFlagChannelActive(false);
        }
    }

    // ---- reset all systematics ---- //
    // loop over all systematics
    for (int isystematic = 0; isystematic < nsystematics; ++isystematic) {
        // get systematic
        BCMTFSystematic* systematic = fMTF->GetSystematic(isystematic);

        // switch off systematic
        if (flag_systematic[isystematic])
            systematic->SetFlagSystematicActive(flag_systematic[isystematic]);
    }

    // ---- reset all channels ---- //

    // loop over all channels
    for (int ichannel = 0; ichannel < nchannels; ++ichannel) {
        // get channel
        BCMTFChannel* channel = fMTF->GetChannel(ichannel);

        // switch channel on/off
        channel->SetFlagChannelActive(flag_channel[ichannel]);
    }

    // ---- workaround: reset MCMC ---- //
    fMTF->ResetResults();

    // ---- print everything ---- //
    TCanvas* c1 = new TCanvas();
    c1->cd();

    // draw first one
    BCMTFComparisonTool* ct =  ctc.at(0);
    ct->DrawOverview();
    //   c1->Print((std::string(filename)+std::string("(")).c_str());
    c1->Print((std::string("overview.pdf") + std::string("(")).c_str());

    // loop over all parameters
    for (int i = 1; i < nparameters - 1; ++ i) {
        // create new comparison tool
        BCMTFComparisonTool* ct = ctc.at(i);

        ct->DrawOverview();
        c1->Print("overview.pdf");
    }

    // draw last one
    ct =  ctc.at(nparameters - 1);
    ct->DrawOverview();
    c1->Print((std::string("overview.pdf") + std::string(")")).c_str());

    // ---- print everything (mcmc) ---- //
    if (flag_mcmc) {
        TCanvas* c2 = new TCanvas();
        c2->cd();

        // draw first one
        BCMTFComparisonTool* ct_mcmc =  ctc_mcmc.at(0);
        ct_mcmc->DrawOverview();
        //   c2->Print((std::string(filename)+std::string("(")).c_str());
        c2->Print((std::string("overview_mcmc.pdf") + std::string("(")).c_str());

        // loop over all parameters
        for (int i = 1; i < nparameters - 1; ++ i) {
            // create new comparison tool
            BCMTFComparisonTool* ct_mcmc = ctc_mcmc.at(i);

            ct_mcmc->DrawOverview();
            c2->Print("overview_mcmc.pdf");
        }

        // draw last one
        ct_mcmc =  ctc_mcmc.at(nparameters - 1);
        ct_mcmc->DrawOverview();
        c2->Print((std::string("overview_mcmc.pdf") + std::string(")")).c_str());
        delete c2;
    }

    // ---- free memory ---- //
    for (int i = 0; i < nparameters; ++i) {
        BCMTFComparisonTool* ct = ctc[i];
        BCMTFComparisonTool* ct_mcmc = ctc_mcmc[i];
        delete ct;
        delete ct_mcmc;
    }
    ctc.clear();
    ctc_mcmc.clear();

    delete c1;

    // ---- change directory ---- //

    cd("../");

    BCLog::OutSummary("Single channel analysis ran successfully");
}

// ---------------------------------------------------------
void BCMTFAnalysisFacility::PerformSingleSystematicAnalyses(const std::string& dirname, const std::string& options)
{
    BCLog::OutSummary("Running single channel systematic analysis in directory \'" + dirname + "\'.");

    // ---- create new directory ---- //

    mkdir(dirname.data(), 0777);
    cd(dirname);

    // ---- check options ---- //

    bool flag_mcmc = true;

    if (options.find("mcmc") < options.size())
        flag_mcmc = true;

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

    // get number of systematics
    int nsystematics = fMTF->GetNSystematics();

    // container of flags for channels
    std::vector<bool> flag_channel(nchannels);

    // container of flags for systematics
    std::vector<bool> flag_systematic(nsystematics);

    // create new container of comparison tools
    std::vector<BCMTFComparisonTool*> ctc;

    // get number of parameters
    int nparameters = fMTF->GetNParameters();

    // ---- add one comparison tool for each systematic ---- //

    // loop over all parameters
    for (int i = 0; i < nparameters; ++ i) {
        // create new comparison tool
        BCMTFComparisonTool* ct = new BCMTFComparisonTool(fMTF->GetParameter(i).GetName());

        // add comparison tool to container
        ctc.push_back(ct);
    }

    // ---- switch on all systematics ---- //

    for (int isystematic = 0; isystematic < nsystematics; ++ isystematic) {
        // get systematic
        BCMTFSystematic* systematic = fMTF->GetSystematic(isystematic);

        // remember old setting
        flag_systematic[isystematic] = systematic->GetFlagSystematicActive();

        // switch on
        systematic->SetFlagSystematicActive(true);
    }

    if (flag_mcmc) {
        // work-around: force initialization
        fMTF->ResetResults();

        // run mcmc
        fMTF->MarginalizeAll();

        // find mode
        fMTF->FindMode(fMTF->GetBestFitParameters());
    } else {
        // find mode
        fMTF->FindMode();
    }

    // print results
    if (flag_mcmc)
        fMTF->PrintAllMarginalized("marginalized_all.pdf");
    fMTF->PrintSummary();

    // loop over all parameters
    for (int i = 0; i < nparameters; ++ i) {
        // get comparison tool
        BCMTFComparisonTool* ct = ctc.at(i);

        ct->AddContribution("all systematics",
                            fMTF->GetBestFitParameters().at(i),
                            fMTF->GetBestFitParameterErrors().at(i));
    }

    // ---- switch off all systematics ---- //

    // loop over all systematics
    for (int isystematic = 0; isystematic < nsystematics; ++isystematic) {
        // get systematic
        BCMTFSystematic* systematic = fMTF->GetSystematic(isystematic);

        // switch off
        systematic->SetFlagSystematicActive(false);
    }

    // ---- perform analysis with all systematics separately ---- //

    // loop over all channels
    for (int isystematic = 0; isystematic < nsystematics; ++isystematic) {
        // get systematic
        BCMTFSystematic* systematic = fMTF->GetSystematic(isystematic);

        // switch on systematic
        systematic->SetFlagSystematicActive(true);

        // perform analysis
        if (flag_mcmc) {
            // work-around: force initialization
            fMTF->ResetResults();

            // run mcmc
            fMTF->MarginalizeAll();

            // find mode
            fMTF->FindMode(fMTF->GetBestFitParameters());
        } else {
            // find mode
            fMTF->FindMode();
        }

        // print results
        if (flag_mcmc)
            fMTF->PrintAllMarginalized(Form("marginalized_systematic_%i.pdf", isystematic));
        BCLog::OutSummary(Form("Summary from fitting systematic %i", isystematic));
        fMTF->PrintFitSummary();

        // ---- update comparison tools ---- //

        // loop over all parameters
        for (int i = 0; i < nparameters; ++ i) {
            // get comparison tool
            BCMTFComparisonTool* ct = ctc.at(i);

            ct->AddContribution(systematic->GetName(),
                                fMTF->GetBestFitParameters().at(i),
                                fMTF->GetBestFitParameterErrors().at(i));
        }

        // switch off systematic
        systematic->SetFlagSystematicActive(false);
    }

    // ---- analysis without any systematic uncertainty ---- //

    if (flag_mcmc) {
        // work-around: force initialization
        fMTF->ResetResults();

        // run mcmc
        fMTF->MarginalizeAll();

        // find mode
        fMTF->FindMode(fMTF->GetBestFitParameters());
    } else {
        // find mode
        fMTF->FindMode();
    }

    // print results
    if (flag_mcmc)
        fMTF->PrintAllMarginalized("marginalized_none.pdf");
    BCLog::OutSummary("Summary without systematics");
    fMTF->PrintFitSummary();

    // loop over all parameters
    for (int i = 0; i < nparameters; ++ i) {
        // get comparison tool
        BCMTFComparisonTool* ct = ctc.at(i);

        ct->AddContribution("no systematics",
                            fMTF->GetBestFitParameters().at(i),
                            fMTF->GetBestFitParameterErrors().at(i));
    }



    // ---- reset all systematics ---- //

    // loop over all systematics
    for (int isystematic = 0; isystematic < nsystematics; ++isystematic) {
        // get systematic
        BCMTFSystematic* systematic = fMTF->GetSystematic(isystematic);

        // switch off systematic
        if (flag_systematic[isystematic])
            systematic->SetFlagSystematicActive(flag_systematic[isystematic]);
    }

    // ---- workaround: reset MCMC ---- //
    fMTF->ResetResults();

    // ---- print everything ---- //
    TCanvas* c1 = new TCanvas();
    c1->cd();

    // draw first one
    BCMTFComparisonTool* ct =  ctc.at(0);
    ct->DrawOverview();
    //   c1->Print((std::string(filename)+std::string("(")).c_str());
    c1->Print((std::string("overview.pdf") + std::string("(")).c_str());

    // loop over all parameters
    for (int i = 1; i < nparameters - 1; ++i) {
        // create new comparison tool
        BCMTFComparisonTool* ct = ctc.at(i);

        ct->DrawOverview();
        c1->Print("overview.pdf");
    }

    // draw last one
    ct =  ctc.at(nparameters - 1);
    ct->DrawOverview();
    c1->Print((std::string("overview.pdf") + std::string(")")).c_str());

    // ---- free memory ---- //
    for (int i = 0; i < nparameters; ++i) {
        BCMTFComparisonTool* ct = ctc[i];
        delete ct;
    }
    ctc.clear();

    delete c1;

    // ---- change directory ---- //

    cd("../");

    BCLog::OutSummary("Single channel analysis ran successfully");
}

// ---------------------------------------------------------
void BCMTFAnalysisFacility::PerformCalibrationAnalysis(const std::string& dirname, const std::vector<double>& default_parameters, int index, const std::vector<double>& parametervalues, int nensembles)
{
    BCLog::OutSummary("Running calibration analysis in directory \'" + dirname + "\'.");

    // ---- create new directory ---- //

    mkdir(dirname.data(), 0777);
    cd(dirname);

    // ---- loop over parameter values and perform analysis  ---- //

    int nvalues = int(parametervalues.size());
    for (int ivalue = 0; ivalue < nvalues; ++ivalue) {

        // open file
        TFile* file = TFile::Open(Form("ensemble_%i.root", ivalue), "RECREATE");
        file->cd();

        // set parameters
        std::vector<double> parameters = default_parameters;
        parameters[index] = parametervalues.at(ivalue);

        // create ensemble
        TTree* tree = PerformEnsembleTest(parameters, nensembles);

        // write tree
        tree->Write();

        // close file
        file->Close();

        // free memory
        delete file;
    }

    // ---- change directory ---- //
    cd("../");

    BCLog::OutSummary("Calibration analysis ran successfully");
}

// ---------------------------------------------------------