BCMTFChannel.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 <iostream>

#include <TCanvas.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TLatex.h>
#include <TMath.h>

#include "BCMTFTemplate.h"
#include "BCMTFSystematicVariation.h"

#include "BCMTFChannel.h"

// ---------------------------------------------------------
BCMTFChannel::BCMTFChannel(const std::string& name)
    : fData(0)
    , fFlagChannelActive(true)
    , fHistUncertaintyBandExpectation(0)
    , fHistUncertaintyBandPoisson(0)
{
    SetName(name);
}

// ---------------------------------------------------------
BCMTFChannel::~BCMTFChannel()
{
    delete fData;

    for (unsigned int i = 0; i < fTemplateContainer.size(); ++i)
        delete (fTemplateContainer.at(i));

    for (unsigned int i = 0; i < fSystematicVariationContainer.size(); ++i)
        delete (fSystematicVariationContainer.at(i));

    delete fHistUncertaintyBandExpectation;
    delete fHistUncertaintyBandPoisson;
}

// ---------------------------------------------------------
void BCMTFChannel::SetHistUncertaintyBandExpectation(TH2D* hist)
{
    if (fHistUncertaintyBandExpectation) {
        delete fHistUncertaintyBandExpectation;
    }
    fHistUncertaintyBandExpectation = hist;
}

// ---------------------------------------------------------
void BCMTFChannel::PrintTemplates(const std::string& filename)
{
    std::string newFilename(filename);
    // check if file extension is pdf
    if ( (newFilename.find_last_of(".") != std::string::npos) &&
            (newFilename.substr(newFilename.find_last_of(".") + 1) == "pdf") ) {
        ; // it's a PDF file

    } else if ( (newFilename.find_last_of(".") != std::string::npos) &&
                (newFilename.substr(newFilename.find_last_of(".") + 1) == "ps") ) {
        ; // it's a PS file
    } else {
        ; // make it a PDF file
        newFilename += ".pdf";
    }

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

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

    int first_hist = -1;
    int last_hist  = 0;

    // calculate first and last existing histogram
    for (int i = 0; i < ntemplates; ++i) {
        // get histogram
        TH1D* temphist = GetTemplate(i)->GetHistogram();

        if (first_hist < 0 && temphist)
            first_hist = i;

        if (temphist)
            last_hist = i;
    }

    // Don't print if there are no histograms
    if ( first_hist < 0 ) {
        // free memory
        delete c1;

        return;
    }

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

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

        TLatex* l = new TLatex();

        // draw
        if (temphist) {
            temphist->Draw();
            l->DrawTextNDC(0.2, 0.9, (fName + "-" + GetTemplate(i)->GetProcessName()).data());
        }

        // print
        if (i == first_hist && (first_hist != last_hist))
            c1->Print((newFilename + "(").data());
        else if (i == last_hist && (first_hist != last_hist))
            c1->Print((newFilename + ")").data());
        else {
            if (temphist)
                c1->Print(newFilename.data());
        }

        // free memory
        delete l;
    }

    // free memory
    delete c1;
}

// ---------------------------------------------------------
void BCMTFChannel::PrintTemplate(int index, const std::string& filename)
{
    // create new canvas
    TCanvas* c1 = new TCanvas();
    c1->cd();

    // get number of systematics
    unsigned int nsystematics = fSystematicVariationContainer.size();

    // draw template
    GetTemplate(index)->GetHistogram()->Draw();

    // print to file
    if (nsystematics == 0)
        c1->Print(filename.data());
    else
        c1->Print( (filename + "(").data() );

    // loop over systematics
    for (unsigned int isystematic = 0; isystematic < nsystematics; ++isystematic) {
        c1->cd();

        // check that histogram exists
        if (!GetSystematicVariation(isystematic)->GetHistogramUp(index) ||
                !GetSystematicVariation(isystematic)->GetHistogramDown(index))
            continue;

        // get histogram
        TH1D hist = TH1D(*GetTemplate(index)->GetHistogram());
        TH1D hist_up(hist);
        TH1D hist_down(hist);

        // set style
        hist.SetFillStyle(0);
        hist_up.SetFillStyle(0);
        hist_up.SetLineStyle(2);
        hist_down.SetFillStyle(0);
        hist_down.SetLineStyle(3);

        // get efficiency
        double efficiency = GetTemplate(index)->GetEfficiency();

        // scale histogram
        hist.Scale(efficiency / hist.Integral());

        // loop over all bins
        for (int i = 1; i <= hist.GetNbinsX(); ++i) {
            hist.SetBinContent(i, GetTemplate(index)->GetHistogram()->GetBinContent(i));
            hist_up.SetBinContent(i, hist.GetBinContent(i) * (1.0 + GetSystematicVariation(isystematic)->GetHistogramUp(index)->GetBinContent(i)));
            hist_down.SetBinContent(i, hist.GetBinContent(i) * (1.0 - GetSystematicVariation(isystematic)->GetHistogramDown(index)->GetBinContent(i)));
        }

        // draw histogram
        hist_up.Draw("HIST");
        hist.Draw("HISTSAME");
        hist_down.Draw("HISTSAME");

        if (isystematic < nsystematics - 1)
            c1->Print(filename.data());
        else
            c1->Print( (filename + ")").data() );
    }

    // free memory
    delete c1;
}

// ---------------------------------------------------------
void BCMTFChannel::PrintHistUncertaintyBandExpectation(const std::string& filename)
{
    // create new canvas
    TCanvas* c1 = new TCanvas();
    c1->cd();

    // draw histogram
    fHistUncertaintyBandExpectation->Draw("COLZ");
    c1->Draw();

    // print
    c1->Print(filename.data());

    // free memory
    delete c1;
}

// ---------------------------------------------------------
void BCMTFChannel::CalculateHistUncertaintyBandPoisson()
{
    // calculate histogram
    int nbinsy_exp = fHistUncertaintyBandExpectation->GetNbinsY();
    int nbinsx_poisson = fHistUncertaintyBandPoisson->GetNbinsX();
    int nbinsy_poisson = fHistUncertaintyBandPoisson->GetNbinsY();

    // loop over x-axis of observation
    for (int ix = 1; ix <= nbinsx_poisson; ++ix) {
        double sum_w = 0;
        // loop over y-axis of expectation and calculate sum of weights
        for (int iy = 1; iy <= nbinsy_exp; ++iy) {
            double w = fHistUncertaintyBandExpectation->GetBinContent(ix, iy);
            sum_w += w;
        }
        // loop over y-axis of expectation
        for (int iy = 1; iy <= nbinsy_exp; ++iy) {
            double w = fHistUncertaintyBandExpectation->GetBinContent(ix, iy) / sum_w;
            double expectation = fHistUncertaintyBandExpectation->GetYaxis()->GetBinCenter(iy);
            // loop over y-axis of observation
            for (int jbin = 1; jbin <= nbinsy_poisson; ++jbin) {
                double p = TMath::Poisson(double(jbin - 1), expectation);
                double bincontent = 0;
                if (iy > 1)
                    bincontent = fHistUncertaintyBandPoisson->GetBinContent(ix, jbin);
                fHistUncertaintyBandPoisson->SetBinContent(ix, jbin, bincontent + p * w);
            }
        }
    }

}

// ---------------------------------------------------------
TH1D* BCMTFChannel::CalculateUncertaintyBandPoisson(double minimum, double maximum, int color)
{
    TH1D* hist;
    {
        BCAux::RootSideEffectGuard g;
        hist = new TH1D(*(fData->GetHistogram()));
    }
    hist->SetMarkerSize(0);
    hist->SetFillColor(color);
    hist->SetFillStyle(1001);

    int nbinsx_poisson = fHistUncertaintyBandPoisson->GetNbinsX();
    int nbinsy_poisson = fHistUncertaintyBandPoisson->GetNbinsY();

    // loop over x-axis of observation
    for (int ix = 1; ix <= nbinsx_poisson; ++ix) {
        double sum_p = 0;  // sum of all probabilities inside the interval
        int limit_min = 0;
        int limit_max = nbinsx_poisson - 1;

        // loop over y-axis of observation
        for (int jbin = 1; jbin <= nbinsy_poisson; ++jbin) {
            double p = fHistUncertaintyBandPoisson->GetBinContent(ix, jbin);
            sum_p += p;
            if (sum_p < minimum)
                limit_min = jbin;
            if (sum_p > maximum && (sum_p - p) < maximum )
                limit_max = jbin - 1;
        }
        //    hist->SetBinContent(ix, 0.5*double(limit_min+limit_max));
        //    hist->SetBinError(ix, 0.5*double(limit_max-limit_min));
        double ylimit_min = fHistUncertaintyBandPoisson->GetYaxis()->GetBinCenter(limit_min);
        double ylimit_max = fHistUncertaintyBandPoisson->GetYaxis()->GetBinCenter(limit_max);
        hist->SetBinContent(ix, 0.5 * double(ylimit_min + ylimit_max));
        hist->SetBinError(ix, 0.5 * double(ylimit_max - ylimit_min));
    }

    return hist;
}

// ---------------------------------------------------------
void BCMTFChannel::PrintHistCumulativeUncertaintyBandPoisson(const std::string& filename)
{
    // create new canvas
    TCanvas* c1 = new TCanvas();
    c1->cd();

    // calculate error band
    this->CalculateHistUncertaintyBandPoisson();

    TH2D hist(*fHistUncertaintyBandPoisson);

    int nbinsx_poisson = hist.GetNbinsX();
    int nbinsy_poisson = hist.GetNbinsY();

    // loop over x-axis of observation
    for (int ix = 1; ix <= nbinsx_poisson; ++ix) {
        double sum_p = 0;  // sum of all probabilities inside the interval

        // loop over y-axis of observation
        for (int jbin = 1; jbin <= nbinsy_poisson; ++jbin) {
            double p = hist.GetBinContent(ix, jbin);
            sum_p += p;
            hist.SetBinContent(ix, jbin, sum_p);
        }
    }

    // draw histogram
    hist.Draw("COLZ");
    c1->Draw();

    // print
    c1->Print(filename.data());

    // free memory
    delete c1;
}

// ---------------------------------------------------------
void BCMTFChannel::PrintHistUncertaintyBandPoisson(const std::string& filename, const std::string& options)
{
    // create new canvas
    TCanvas* c1 = new TCanvas();
    c1->cd();

    // calculate error band
    this->CalculateHistUncertaintyBandPoisson();

    // draw histogram
    fHistUncertaintyBandPoisson->Draw(options.c_str());
    c1->Draw();

    // print
    c1->Print(filename.data());

    // free memory
    delete c1;
}

// ---------------------------------------------------------
void BCMTFChannel::PrintUncertaintyBandPoisson(const std::string& filename, double minimum, double maximum, int color)
{
    // create new canvas
    TCanvas* c1 = new TCanvas();
    c1->cd();

    // calculate error band
    TH1D* hist = CalculateUncertaintyBandPoisson(minimum, maximum, color);

    // draw histogram
    hist->Draw("E2");
    c1->Draw();

    // print
    c1->Print(filename.data());

    // free memory
    delete c1;
}

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