Lau1DHistPdf.cc

Go to the documentation of this file.

// Copyright University of Warwick 2006 - 2013.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

// Authors:
// Thomas Latham
// John Back
// Paul Harrison

#include <cstdlib>
#include <iostream>
#include <vector>

#include "TAxis.h"
#include "TH1.h"
#include "TRandom.h"
#include "TSystem.h"

#include "Lau1DHistPdf.hh"
#include "LauRandom.hh"

class LauParameter;

ClassImp(Lau1DHistPdf)


Lau1DHistPdf::Lau1DHistPdf(const TString& theVarName, const TH1* hist, Double_t minAbscissa, Double_t maxAbscissa,
                Bool_t useInterpolation, Bool_t fluctuateBins) :
        LauAbsPdf(theVarName, std::vector<LauParameter*>(), minAbscissa, maxAbscissa),
        hist_(hist ? dynamic_cast<TH1*>(hist->Clone()) : 0),
        useInterpolation_(useInterpolation),
        fluctuateBins_(fluctuateBins),
        nBins_(0),
        axisMin_(0.0),
        axisMax_(0.0),
        axisRange_(0.0)
{
        // Constructor

        // Set the directory for the histogram
        hist_->SetDirectory(0);

        // Save various attributes of the histogram 
        nBins_ = hist_->GetNbinsX();
        TAxis* xAxis = hist_->GetXaxis();
        axisMin_ = xAxis->GetXmin();
        axisMax_ = xAxis->GetXmax();
        axisRange_ = axisMax_ - axisMin_;

        // Check that axis range corresponds to range of abscissa
        if (TMath::Abs(this->getMinAbscissa() - axisMin_)>1e-6) {
                std::cerr << "ERROR in Lau1DHistPdf::Lau1DHistPdf : Histogram axis minimum: " << axisMin_ << 
                        " does not correspond to abscissa minimum: " << this->getMinAbscissa() << "." << std::endl;
                gSystem->Exit(EXIT_FAILURE);
        }
        if (TMath::Abs(this->getMaxAbscissa() - axisMax_)>1e-6) {
                std::cerr << "ERROR in Lau1DHistPdf::Lau1DHistPdf : Histogram axis maximum: " << axisMax_ << 
                        " does not correspond to abscissa maximum: " << this->getMaxAbscissa() << "." << std::endl;
                gSystem->Exit(EXIT_FAILURE);
        }

        // If the bins are to be fluctuated then do so now before
        // calculating anything that depends on the bin content.
        if (fluctuateBins) {
                this->doBinFluctuation();
        }

        // Calculate the PDF normalisation.
        this->calcNorm();

        // And check it is OK.
        this->checkNormalisation();
}

Lau1DHistPdf::~Lau1DHistPdf()
{
        // Destructor
        delete hist_; hist_ = 0;
}

Lau1DHistPdf::Lau1DHistPdf(const Lau1DHistPdf& other) : LauAbsPdf(other.varName(), other.getParameters(), other.getMinAbscissa(), other.getMaxAbscissa())
{
        hist_ = other.hist_ ? dynamic_cast<TH1*>(other.hist_->Clone()) : 0;
        useInterpolation_ = other.useInterpolation_;
        fluctuateBins_ = other.fluctuateBins_;
        this->setRandomFun(other.getRandomFun());
        this->calcNorm();
}

void Lau1DHistPdf::calcPDFHeight( const LauKinematics* /*kinematics*/ )
{
        if (this->heightUpToDate()) {
                return;
        }

        // Get the maximum height of the histogram
        Int_t maxBin = hist_->GetMaximumBin();
        Double_t height = hist_->GetBinContent(maxBin);
        this->setMaxHeight(height);
}

void Lau1DHistPdf::calcNorm()
{
        // Calculate the histogram normalisation.

        // Loop over the range to get the total area.
        // Just sum the contributions up using 1e-3 increments of the range.
        // Multiply the end result by dx.

        Double_t dx(1e-3*axisRange_);
        Double_t area(0.0);

        Double_t x(axisMin_ + dx/2.0);
        while (x > axisMin_ && x < axisMax_) {
                area += this->interpolate(x);
                x += dx;
        }

        Double_t norm = area*dx;
        this->setNorm(norm);
}

void Lau1DHistPdf::checkNormalisation()
{
        Double_t dx(1e-3*axisRange_);
        Double_t area(0.0);
        Double_t areaNoNorm(0.0);

        Double_t x(axisMin_ + dx/2.0);
        while (x > axisMin_ && x < axisMax_) {
                area += this->interpolateNorm(x);
                useInterpolation_ = kFALSE;
                areaNoNorm += this->interpolate(x);
                useInterpolation_ = kTRUE;
                x += dx;
        }
        Double_t norm = area*dx;

        std::cout << "INFO in Lau1DHistPdf::checkNormalisation : Area = " << area << ", dx = " << dx << std::endl;
        std::cout << "                                         : Area with no norm = " << areaNoNorm << "*dx = " << areaNoNorm*dx << std::endl;
        std::cout << "                                         : The total area of the normalised histogram PDF is " << norm << std::endl;
}

Double_t Lau1DHistPdf::getBinHistValue(Int_t bin) const
{
        // Check that bin is in range [1 , nBins_]
        if ((bin < 1) ||  (bin > nBins_)) {
                return 0.0;
        }
        Double_t value = static_cast<Double_t>(hist_->GetBinContent(bin));
        // protect against negative values
        if ( value < 0.0 ) {
                std::cerr << "WARNING in Lau1DHistPdf::getBinHistValue : Negative bin content set to zero!" << std::endl;
                value = 0.0;
        }
        return value;
}

Double_t Lau1DHistPdf::interpolateNorm(Double_t x) const
{
        // Get the normalised interpolated value.
        Double_t value = this->interpolate(x);
        Double_t norm = this->getNorm();
        return value/norm;
}

void Lau1DHistPdf::calcLikelihoodInfo(const LauAbscissas& abscissas)
{
        // Check that the given abscissa is within the allowed range
        if (!this->checkRange(abscissas)) {
                gSystem->Exit(EXIT_FAILURE);
        }

        // Get our abscissa
        Double_t abscissa = abscissas[0];

        // Calculate the interpolated value
        Double_t value = this->interpolate(abscissa);
        this->setUnNormPDFVal(value);
}

Double_t Lau1DHistPdf::interpolate(Double_t x) const
{
        // This function returns the interpolated value of the histogram function
        // for the given value of x by finding the adjacent bins and extrapolating
        // using weights based on the inverse distance of the point from the adajcent
        // bin centres.

        // Find the histogram bin
        Int_t bin = hist_->FindBin(x);

        // Ask whether we want to do the interpolation, since this method is
        // not reliable for low statistics histograms.
        if (useInterpolation_ == kFALSE) {
                return this->getBinHistValue(bin);
        }

        // Find the bin centres (actual co-ordinate positions, not histogram indices)
        Double_t cbinx = hist_->GetBinCenter(bin);

        // Find the adjacent bins
        Double_t deltax = x - cbinx;

        Int_t bin_adj(0);
        if (deltax > 0.0) {
                bin_adj = bin + 1;
        } else {
                bin_adj = bin - 1;
        }


        Bool_t isBoundary(kFALSE);
        if ( bin_adj > nBins_ || bin_adj < 1 ) {
                isBoundary = kTRUE;
        }

        // At the edges, do no interpolation, use entry in bin.
        if (isBoundary == kTRUE) {
                return this->getBinHistValue(bin);
        }

        // Linear interpolation using inverse distance as weights.
        // Find the adjacent bin centre
        Double_t cbinx_adj = hist_->GetBinCenter(bin_adj);
        Double_t deltax_adj = cbinx_adj - x;

        Double_t dx0 = TMath::Abs(deltax);
        Double_t dx1 = TMath::Abs(deltax_adj);

        Double_t denom = dx0 + dx1;

        Double_t value0 = this->getBinHistValue(bin);
        Double_t value1 = this->getBinHistValue(bin_adj);

        Double_t value = value0*dx1 + value1*dx0;
        value /= denom;

        return value;
}

void Lau1DHistPdf::doBinFluctuation()
{
        TRandom* random = LauRandom::zeroSeedRandom();
        for (Int_t bin(0); bin<nBins_; bin++) {
                Double_t currentContent = hist_->GetBinContent(bin+1);
                Double_t currentError   = hist_->GetBinError(bin+1);
                Double_t newContent = random->Gaus(currentContent,currentError);
                if (newContent<0.0) {
                        hist_->SetBinContent(bin+1,0.0);
                } else {
                        hist_->SetBinContent(bin+1,newContent);
                }
        }
}