Lau2DAbsHistDP.cc

Go to the documentation of this file.

// Copyright University of Warwick 2004 - 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 <iostream>

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

#include "Lau2DAbsHistDP.hh"
#include "LauBifurcatedGaussPdf.hh"
#include "LauDaughters.hh"
#include "LauKinematics.hh"
#include "LauRandom.hh"

ClassImp(Lau2DAbsHistDP)


Lau2DAbsHistDP::Lau2DAbsHistDP(const LauDaughters* daughters, Bool_t useUpperHalfOnly, Bool_t squareDP) :
        kinematics_( (daughters!=0) ? daughters->getKinematics() : 0 ),
        upperHalf_(useUpperHalfOnly),
        squareDP_(squareDP)
{
        if ( squareDP && ! daughters->squareDP() ) {
                // The histogram provided is defined in the square DP but the
                // kinematics object has calculation of the square DP
                // co-ordinates disabled, so need to enable it,
                // which requires a bit of a unpleasant const_cast
                std::cerr << "WARNING in Lau2DAbsHistDP constructor : forcing kinematics to calculate the required square DP co-ordinates" << std::endl;
                LauKinematics* kine = const_cast<LauKinematics*>( kinematics_ );
                kine->squareDP(kTRUE);
        }
}

Lau2DAbsHistDP::~Lau2DAbsHistDP()
{
}

void Lau2DAbsHistDP::doBinFluctuation(TH2* hist)
{
        TRandom* random = LauRandom::randomFun();

        Int_t nBinsX = static_cast<Int_t>(hist->GetNbinsX());
        Int_t nBinsY = static_cast<Int_t>(hist->GetNbinsY());

        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);
                        }
                }
        }
}

void Lau2DAbsHistDP::doBinFluctuation(TH2* hist, const TH2* errorHi, const TH2* errorLo)
{
        LauParameter* mean = new LauParameter("mean",   0.5, 0.0, 1.0, kFALSE);
        LauParameter* sigL = new LauParameter("sigmaL", 0.5, 0.0, 1.0, kFALSE);
        LauParameter* sigR = new LauParameter("sigmaR", 0.5, 0.0, 1.0, kFALSE);

        std::vector<LauAbsRValue*> pars(3);
        pars[0] = mean;
        pars[1] = sigL;
        pars[2] = sigR;

        const TString varName("tmp");

        Int_t nBinsX = static_cast<Int_t>(hist->GetNbinsX());
        Int_t nBinsY = static_cast<Int_t>(hist->GetNbinsY());

        LauFitData genData;

        for (Int_t i(0); i<nBinsX; ++i) {
                for (Int_t j(0); j<nBinsY; ++j) {
                        const Double_t currentContent = hist->GetBinContent(i+1,j+1);
                        const Double_t currentErrorLo = errorLo->GetBinContent(i+1,j+1);
                        const Double_t currentErrorHi = errorHi->GetBinContent(i+1,j+1);

                        mean->value( currentContent );
                        sigL->value( currentErrorLo );
                        sigR->value( currentErrorHi );

                        const Double_t minVal = TMath::Max( 0.0, currentContent-5.0*currentErrorLo );
                        const Double_t maxVal = TMath::Min( 1.0, currentContent+5.0*currentErrorHi );

                        LauBifurcatedGaussPdf bfgaus(varName, pars, minVal, maxVal);
                        bfgaus.heightUpToDate(kFALSE);
                        genData = bfgaus.generate(0);

                        const Double_t newContent = genData[varName];

                        hist->SetBinContent(i+1,j+1,newContent);
                }
        }

        delete pars[0];
        delete pars[1];
        delete pars[2];
        pars.clear();
}

void Lau2DAbsHistDP::raiseOrLowerBins(TH2* hist, const Double_t avEff, const Double_t avEffError)
{
        TRandom* random = LauRandom::randomFun();

        Double_t curAvg = this->computeAverageContents(hist);
        Double_t newAvg = random->Gaus(avEff,avEffError);

        hist->Scale( newAvg / curAvg );
}

Double_t Lau2DAbsHistDP::computeAverageContents(const TH2* hist) const
{
        Double_t totalContent(0.0);
        Int_t binsWithinDPBoundary(0);

        Int_t nBinsX = static_cast<Int_t>(hist->GetNbinsX());
        Int_t nBinsY = static_cast<Int_t>(hist->GetNbinsY());

        // Loop through the bins and include any that have their centre or any
        // of the four corners within the kinematic boundary
        for ( Int_t i(0); i<nBinsX; ++i ) {
                Double_t binXCentre = hist->GetXaxis()->GetBinCenter(i+1);
                Double_t binXLowerEdge = hist->GetXaxis()->GetBinLowEdge(i+1);
                Double_t binXUpperEdge = hist->GetXaxis()->GetBinUpEdge(i+1);

                for ( Int_t j(0); j<nBinsY; ++j ) {
                        Double_t binYCentre = hist->GetYaxis()->GetBinCenter(i+1);
                        Double_t binYLowerEdge = hist->GetYaxis()->GetBinLowEdge(i+1);
                        Double_t binYUpperEdge = hist->GetYaxis()->GetBinUpEdge(i+1);

                        if ( this->withinDPBoundaries( binXCentre, binYCentre ) ||
                             this->withinDPBoundaries( binXLowerEdge, binYLowerEdge ) ||
                             this->withinDPBoundaries( binXUpperEdge, binYUpperEdge ) ||
                             this->withinDPBoundaries( binXLowerEdge, binYUpperEdge ) ||
                             this->withinDPBoundaries( binXUpperEdge, binYLowerEdge ) ) {

                                totalContent += hist->GetBinContent(i+1, j+1);
                                ++binsWithinDPBoundary;
                        }
                }
        }

        return totalContent/binsWithinDPBoundary;
}

Bool_t Lau2DAbsHistDP::withinDPBoundaries(Double_t x, Double_t y) const
{
        return squareDP_ ? kinematics_->withinSqDPLimits(x,y) : kinematics_->withinDPLimits(x,y);
}

void Lau2DAbsHistDP::getUpperHalf(Double_t& x, Double_t& y) const
{
        if ( upperHalf_ == kTRUE ) {
                if ( squareDP_ == kFALSE && x > y ) {
                        Double_t temp = y;
                        y = x;
                        x = temp;
                } else if ( squareDP_ == kTRUE && y > 0.5 ) {
                        y = 1.0 - y;
                }
        }
}