Lau2DHistPdf.cc

Go to the documentation of this file.

// Copyright University of Warwick 2008 - 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 "TH2.h"
#include "TRandom.h"
#include "TSystem.h"

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

class LauParameter;

ClassImp(Lau2DHistPdf)


Lau2DHistPdf::Lau2DHistPdf(const std::vector<TString>& theVarNames, const TH2* hist,
                const LauFitData& minVals, const LauFitData& maxVals,
                Bool_t useInterpolation, Bool_t fluctuateBins) :
        LauAbsPdf(theVarNames, std::vector<LauAbsRValue*>(), minVals, maxVals),
        hist_(hist ? dynamic_cast<TH2*>(hist->Clone()) : 0),
        xProj_(0),
        yProj_(0),
        xVarPdf_(0),
        yVarPdf_(0),
        xName_(""),
        yName_(""),
        nBinsX_(0),
        nBinsY_(0),
        minX_(0.0),
        maxX_(0.0),
        minY_(0.0),
        maxY_(0.0),
        rangeX_(0.0),
        rangeY_(0.0),
        binXWidth_(0.0),
        binYWidth_(0.0),
        invBinXWidth_(0.0),
        invBinYWidth_(0.0),
        useInterpolation_(useInterpolation),
        fluctuateBins_(fluctuateBins)
{
        // Constructor

        // Check that we have two variables
        if ( this->nInputVars() != 2 ) {
                std::cerr << "ERROR in Lau2DHistPdf::Lau2DHistPdf : Have not been provided with exactly two variables." << std::endl;
                gSystem->Exit(EXIT_FAILURE);
        }

        // Set the variable names from the abstract class
        xName_ = this->varNames()[0];
        yName_ = this->varNames()[1];

        // Set the variable limits from the abstract class
        minX_ = this->getMinAbscissa(xName_);
        maxX_ = this->getMaxAbscissa(xName_);
        minY_ = this->getMinAbscissa(yName_);
        maxY_ = this->getMaxAbscissa(yName_);
        rangeX_ = maxX_ - minX_;
        rangeY_ = maxY_ - minY_;

        // Have we got a valid histogram
        if ( hist_ == 0 ) {
                std::cerr << "ERROR in Lau2DHistPdf::Lau2DHistPdf : Histogram pointer is null." << std::endl;
                gSystem->Exit(EXIT_FAILURE);
        }

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

        // Save various attributes of the histogram 
        nBinsX_ = hist_->GetNbinsX();
        nBinsY_ = hist_->GetNbinsY();

        TAxis* xAxis = hist_->GetXaxis();
        Double_t xAxisMin = xAxis->GetXmin();
        Double_t xAxisMax = xAxis->GetXmax();

        TAxis* yAxis = hist_->GetYaxis();
        Double_t yAxisMin = yAxis->GetXmin();
        Double_t yAxisMax = yAxis->GetXmax();

        // Check that axis ranges corresponds to ranges of abscissas
        if (TMath::Abs(minX_ - xAxisMin)>1e-6) {
                std::cerr << "ERROR in Lau2DHistPdf::Lau2DHistPdf : Histogram x-axis minimum: " << xAxisMin <<
                             " does not correspond to " << xName_ << " minimum: " << minX_ << "." << std::endl;
                gSystem->Exit(EXIT_FAILURE);
        }
        if (TMath::Abs(maxX_ - xAxisMax)>1e-6) {
                std::cerr << "ERROR in Lau2DHistPdf::Lau2DHistPdf : Histogram x-axis maximum: " << xAxisMax <<
                        " does not correspond to " << xName_ << " maximum: " << maxX_ << "." << std::endl;
                gSystem->Exit(EXIT_FAILURE);
        }
        if (TMath::Abs(minY_ - yAxisMin)>1e-6) {
                std::cerr << "ERROR in Lau2DHistPdf::Lau2DHistPdf : Histogram y-axis minimum: " << yAxisMin <<
                        " does not correspond to " << yName_ << " minimum: " << minY_ << "." << std::endl;
                gSystem->Exit(EXIT_FAILURE);
        }
        if (TMath::Abs(maxY_ - yAxisMax)>1e-6) {
                std::cerr << "ERROR in Lau2DHistPdf::Lau2DHistPdf : Histogram y-axis maximum: " << yAxisMax <<
                        " does not correspond to " << yName_ << " maximum: " << maxY_ << "." << std::endl;
                gSystem->Exit(EXIT_FAILURE);
        }

        // Calculate the bin widths and inverse bin widths
        binXWidth_ = static_cast<Double_t>(TMath::Abs(rangeX_)/(nBinsX_*1.0));
        binYWidth_ = static_cast<Double_t>(TMath::Abs(rangeY_)/(nBinsY_*1.0));
        if (binXWidth_ > 1e-10) {invBinXWidth_ = 1.0/binXWidth_;}
        if (binYWidth_ > 1e-10) {invBinYWidth_ = 1.0/binYWidth_;}

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

        // Create the projections and 1D PDFs from those
        xProj_ = hist_->ProjectionX();
        yProj_ = hist_->ProjectionY();
        if ( !xProj_ || !yProj_ ) {
                std::cerr << "ERROR in Lau2DHistPdf::Lau2DHistPdf : Can't get X or Y projection from 2D histogram." << std::endl;
                gSystem->Exit(EXIT_FAILURE);
        }
        xVarPdf_ = new Lau1DHistPdf(xName_, xProj_, minX_, maxX_, useInterpolation_, fluctuateBins_);
        yVarPdf_ = new Lau1DHistPdf(yName_, yProj_, minY_, maxY_, useInterpolation_, fluctuateBins_);

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

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

Lau2DHistPdf::~Lau2DHistPdf()
{
        // Destructor
        delete hist_; hist_ = 0;
        delete xProj_; xProj_ = 0;
        delete yProj_; yProj_ = 0;
        delete xVarPdf_; xVarPdf_ = 0;
        delete yVarPdf_; yVarPdf_ = 0;
}

void Lau2DHistPdf::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 Lau2DHistPdf::calcNorm()
{
        // Calculate the histogram normalisation.

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

        Double_t dx(1e-3*rangeX_), dy(1e-3*rangeY_);
        Double_t area(0.0);

        Double_t x(minX_ + dx/2.0);
        while (x > minX_ && x < maxX_) {
                Double_t y(minY_ + dy/2.0);
                while (y > minY_ && y < maxY_) {        
                        area += this->interpolateXY(x,y);
                        y += dy;
                } // y while loop
                x += dx;
        } // x while loop

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

void Lau2DHistPdf::checkNormalisation()
{
        Double_t dx(1e-3*rangeX_), dy(1e-3*rangeY_);
        Double_t area(0.0);
        Double_t areaNoNorm(0.0);

        // Preserve the value of a variable we change temporarily
        Bool_t interpolate = useInterpolation_;

        Double_t x(minX_ + dx/2.0);
        while (x > minX_ && x < maxX_) {
                Double_t y(minY_ + dy/2.0);
                while (y > minY_ && y < maxY_) {        
                        area += this->interpolateXYNorm(x,y);
                        useInterpolation_ = kFALSE;
                        areaNoNorm += this->interpolateXY(x,y);
                        useInterpolation_ = kTRUE;
                        y += dy;
                } // y while loop
                x += dx;
        } // x while loop
        Double_t norm = area*dx*dy;

        useInterpolation_ = interpolate; //Restore old value

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

Double_t Lau2DHistPdf::getBinHistValue(Int_t i, Int_t j) const
{
        // Check that x co-ord is in range [1 , nBinsX_]
        if ((i < 1) ||  (i > nBinsX_)) {
                return 0.0;
        }
        // Check that y co-ord is in range [1 , nBinsY_]
        if ((j < 1) ||  (j > nBinsY_)) {
                return 0.0;
        }
        Double_t value = hist_->GetBinContent(i,j);
        // protect against negative values
        if ( value < 0.0 ) {
                std::cerr << "WARNING in Lau2DHistPdf::getBinHistValue : Negative bin content set to zero!" << std::endl;
                value = 0.0;
        }
        return value;
}

Double_t Lau2DHistPdf::interpolateXYNorm(Double_t x, Double_t y) const
{
        // Get the normalised interpolated value.
        Double_t value = this->interpolateXY(x,y);
        Double_t norm = this->getNorm();
        return value/norm;
}

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

        // Find the 2D histogram bin for x and y
        Int_t i = static_cast<Int_t>((x - minX_)*invBinXWidth_)+1;
        Int_t j = static_cast<Int_t>((y - minY_)*invBinYWidth_)+1;

        if (i > nBinsX_) {i = nBinsX_;}
        if (j > nBinsY_) {j = nBinsY_;}

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

        // Find the bin centres (actual co-ordinate positions, not histogram indices)
        Double_t cbinx = static_cast<Double_t>(i-0.5)*binXWidth_ + minX_;
        Double_t cbiny = static_cast<Double_t>(j-0.5)*binYWidth_ + minY_;

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

        Int_t i_adj(0), j_adj(0);
        if (deltax > 0.0) {
                i_adj = i + 1;
        } else {
                i_adj = i - 1;
        }
        if (deltay > 0.0) {
                j_adj = j + 1;
        } else {
                j_adj = j - 1;
        }

        Bool_t isXBoundary(kFALSE), isYBoundary(kFALSE);

        Double_t value(0.0);

        if (i_adj > nBinsX_ || i_adj < 1) {isYBoundary = kTRUE;}
        if (j_adj > nBinsY_ || j_adj < 1) {isXBoundary = kTRUE;}

        // In the corners, do no interpolation. Use entry in bin (i,j)
        if (isXBoundary == kTRUE && isYBoundary == kTRUE) {

                value = this->getBinHistValue(i,j);

        } else if (isXBoundary == kTRUE && isYBoundary == kFALSE) {

                // Find the adjacent x bin centre
                Double_t cbinx_adj = static_cast<Double_t>(i_adj-0.5)*binXWidth_ + minX_;

                Double_t dx0 = TMath::Abs(x - cbinx);
                Double_t dx1 = TMath::Abs(cbinx_adj - x);
                Double_t inter_denom = dx0 + dx1;

                Double_t value1 = this->getBinHistValue(i,j);
                Double_t value2 = this->getBinHistValue(i_adj,j);

                value = (value1*dx1 + value2*dx0)/inter_denom;

        } else if (isYBoundary == kTRUE && isXBoundary == kFALSE) {

                // Find the adjacent y bin centre
                Double_t cbiny_adj = static_cast<Double_t>(j_adj-0.5)*binYWidth_ + minY_;

                Double_t dy0 = TMath::Abs(y - cbiny);
                Double_t dy1 = TMath::Abs(cbiny_adj - y);
                Double_t inter_denom = dy0 + dy1;

                Double_t value1 = this->getBinHistValue(i,j);
                Double_t value2 = this->getBinHistValue(i,j_adj);

                value = (value1*dy1 + value2*dy0)/inter_denom;

        } else {

                // 2D linear interpolation using inverse distance as weights.
                // Find the adjacent x and y bin centres
                Double_t cbinx_adj = static_cast<Double_t>(i_adj-0.5)*binXWidth_ + minX_;
                Double_t cbiny_adj = static_cast<Double_t>(j_adj-0.5)*binYWidth_ + minY_;

                Double_t dx0 = TMath::Abs(x - cbinx);
                Double_t dx1 = TMath::Abs(cbinx_adj - x);
                Double_t dy0 = TMath::Abs(y - cbiny);
                Double_t dy1 = TMath::Abs(cbiny_adj - y);

                Double_t inter_denom = (dx0 + dx1)*(dy0 + dy1);

                Double_t value1 = this->getBinHistValue(i,j);
                Double_t value2 = this->getBinHistValue(i_adj,j);
                Double_t value3 = this->getBinHistValue(i,j_adj);
                Double_t value4 = this->getBinHistValue(i_adj,j_adj);

                value = value1*dx1*dy1 + value2*dx0*dy1 + value3*dx1*dy0 + value4*dx0*dy0;
                value /= inter_denom;
        }

        return value;
}

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

        // Calculate the interpolated value
        Double_t x = abscissas[0];
        Double_t y = abscissas[1];
        Double_t value = this->interpolateXY(x,y);
        this->setUnNormPDFVal(value);

        // TODO - do we need this?  I think probably not.
        // Have the 1D PDFs calculate their likelihoods as well
        //xVarPdf_->calcLikelihoodInfo(x);
        //yVarPdf_->calcLikelihoodInfo(y);
}

void Lau2DHistPdf::calcLikelihoodInfo(UInt_t iEvt)
{
        // Call the base class method
        this->LauAbsPdf::calcLikelihoodInfo(iEvt);

        // Have the 1D PDFs retrieve their values as well
        xVarPdf_->calcLikelihoodInfo(iEvt);
        yVarPdf_->calcLikelihoodInfo(iEvt);
}

Double_t Lau2DHistPdf::getLikelihood( const TString& theVarName ) const
{
        if ( theVarName == xName_ ) {
                return xVarPdf_->getLikelihood();
        } else if ( theVarName == yName_ ) {
                return yVarPdf_->getLikelihood();
        } else {
                std::cerr << "ERROR in Lau2DHistPdf::getLikelihood : Unrecognised variable name \"" << theVarName << "\", cannot determine likelihood." << std::endl;
                return 0.0;
        }
}

void Lau2DHistPdf::cacheInfo(const LauFitDataTree& inputData)
{
        // Call the base class method
        this->LauAbsPdf::cacheInfo(inputData);

        // Have the 1D PDFs cache their info as well
        xVarPdf_->cacheInfo(inputData);
        yVarPdf_->cacheInfo(inputData);
}

void Lau2DHistPdf::doBinFluctuation()
{
        TRandom* random(LauRandom::randomFun());
        for (Int_t i(0); i<nBinsX_; i++) {
                for (Int_t j(0); j<nBinsY_; j++) {
                        Double_t currentContent = hist_->GetBinContent(i+1,j+1);
                        Double_t currentError   = hist_->GetBinError(i+1,j+1);
                        Double_t newContent = random->Gaus(currentContent,currentError);
                        if (newContent<0.0) {
                                hist_->SetBinContent(i+1,j+1,0.0);
                        } else {
                                hist_->SetBinContent(i+1,j+1,newContent);
                        }
                }
        }
}

LauFitData Lau2DHistPdf::generate(const LauKinematics* kinematics)
{
        this->withinGeneration(kTRUE);

        // Check that the PDF height is up to date
        this->calcPDFHeight( kinematics );

        // Generate the value of the abscissa.
        Bool_t gotAbscissa(kFALSE);
        if (this->getRandomFun() == 0) {
                std::cerr << "ERROR in Lau2DHistPdf::generate : Please set the random number generator for this PDF by using the setRandomFun(TRandom*) function." << std::endl;
                this->withinGeneration(kFALSE);
                gSystem->Exit(EXIT_FAILURE);
        }

        Double_t genPDFVal(0.0);
        LauAbscissas genAbscissas(2);

        Double_t PDFheight = this->getMaxHeight()*(1.0+1e-11);
        while (!gotAbscissa) {

                genAbscissas[0] = this->getRandomFun()->Rndm()*this->getRange(xName_) + this->getMinAbscissa(xName_);
                genAbscissas[1] = this->getRandomFun()->Rndm()*this->getRange(yName_) + this->getMinAbscissa(yName_);

                this->calcLikelihoodInfo(genAbscissas);
                genPDFVal = this->getUnNormLikelihood();

                if (this->getRandomFun()->Rndm() <= genPDFVal/PDFheight) {gotAbscissa = kTRUE;}

                if (genPDFVal > PDFheight) {
                        std::cerr << "WARNING in LauAbsPdf::generate : genPDFVal = " << genPDFVal << " is larger than the specified PDF height " << this->getMaxHeight() << " for (x,y) = (" << genAbscissas[0] << "," << genAbscissas[1] << ")." << std::endl;
                        std::cerr << "                               : Need to reset height to be larger than " << genPDFVal << " by using the setMaxHeight(Double_t) function and re-run the Monte Carlo generation!" << std::endl;
                }
        }

        LauFitData genData;
        genData[ xName_ ] = genAbscissas[0];
        genData[ yName_ ] = genAbscissas[1];

        this->withinGeneration(kFALSE);

        return genData;
}