LauDPDepBifurGaussPdf.cc

Go to the documentation of this file.

// Copyright University of Warwick 2007 - 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 <vector>
using std::cout;
using std::cerr;
using std::endl;
using std::vector;

#include "TMath.h"
#include "TRandom.h"
#include "TSystem.h"

#include "Lau1DHistPdf.hh"
#include "LauConstants.hh"
#include "LauComplex.hh"
#include "LauDaughters.hh"
#include "LauFitDataTree.hh"
#include "LauKinematics.hh"
#include "LauDPDepBifurGaussPdf.hh"

ClassImp(LauDPDepBifurGaussPdf)


LauDPDepBifurGaussPdf::LauDPDepBifurGaussPdf(const TString& theVarName, const vector<LauAbsRValue*>& params,
                Double_t minAbscissa, Double_t maxAbscissa,
                const LauDaughters* daughters,
                const vector<Double_t>& meanCoeffs,
                const vector<Double_t>& sigmaLCoeffs,
                const vector<Double_t>& sigmaRCoeffs,
                DPAxis dpAxis) :
        LauAbsPdf(theVarName, params, minAbscissa, maxAbscissa),
        kinematics_(daughters ? daughters->getKinematics() : 0),
        mean_(0),
        sigmaL_(0),
        sigmaR_(0),
        meanVal_(0.0),
        sigmaLVal_(0.0),
        sigmaRVal_(0.0),
        meanCoeffs_(meanCoeffs),
        sigmaLCoeffs_(sigmaLCoeffs),
        sigmaRCoeffs_(sigmaRCoeffs),
        dpAxis_(dpAxis),
        scaleMethod_(poly)
{
        // Check we have a valid kinematics object
        if ( ! kinematics_ ) {
                cerr<<"ERROR in LauDPDepBifurGaussPdf::LauDPDepBifurGaussPdf : Have not been provided with a valid DP kinematics object."<<endl;
                gSystem->Exit(EXIT_FAILURE);
        }

        // The parameters are:
        // - the mean, the sigmaL and the sigmaR of the Bifurcated Gaussian
        //
        // The next two arguments specify the range in which the PDF is defined,
        // and the PDF will be normalised w.r.t. these limits.
        //
        // The final three argument define whether the BF Gaussian parameters should be scaled or not

        mean_ = this->findParameter("mean");
        sigmaL_ = this->findParameter("sigmaL");
        sigmaR_ = this->findParameter("sigmaR");

        if ((this->nParameters() != 3) || (mean_ == 0) || (sigmaL_ == 0) || (sigmaR_ == 0) ) {
                cerr<<"ERROR in LauDPDepBifurGaussPdf constructor: LauDPDepBifurGaussPdf requires 3 parameters:"
                        <<" \"mean\", \"sigmaL\" and \"sigmaR\"."<<endl;
                gSystem->Exit(EXIT_FAILURE);
        }

        // Cache the normalisation factor
        this->calcNorm();
}

LauDPDepBifurGaussPdf::~LauDPDepBifurGaussPdf() 
{
        // Destructor
}

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

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

        // Get the up to date parameter values
        meanVal_ = mean_->unblindValue();
        sigmaLVal_ = sigmaL_->unblindValue();
        sigmaRVal_ = sigmaR_->unblindValue();
        
        // Find out the DP position
        Double_t dpPos(0.0);
        UInt_t nVars = this->nInputVars();
        if ( abscissas.size() == nVars+2 ) {
                Double_t m13Sq = abscissas[nVars];
                Double_t m23Sq = abscissas[nVars+1];

                if ( dpAxis_ == M12 ) {
                        dpPos = kinematics_->calcThirdMassSq(m13Sq,m23Sq);
                } else if ( dpAxis_ == M13 ) {
                        dpPos = m13Sq;
                } else if ( dpAxis_ == M23 ) {
                        dpPos = m23Sq;
                } else if ( dpAxis_ == MMIN ) {
                        dpPos = TMath::Min( m13Sq, m23Sq );
                } else if ( dpAxis_ == MMAX ) {
                        dpPos = TMath::Max( m13Sq, m23Sq );
                } else {
                        dpPos = kinematics_->distanceFromDPCentre(m13Sq,m23Sq);
                }
        }

        // Scale the gaussian parameters by the dpPos (if appropriate)
        ScaleMethod scale = this->scaleMethod();
        if ( scale==poly ) {
                this->scalePars_poly(dpPos);
        } else if ( scale==polyNegPower ) {
                this->scalePars_polyNegPower(dpPos);
        } else {
                cerr<<"Scaling method unknown!  Methods available: <poly> and <polyNegPower>."<<endl;
                gSystem->Exit(EXIT_FAILURE);
        }

        // Evaluate the Birfucated Gaussian PDF value
        Double_t arg = abscissa - meanVal_;
        Double_t coef(0.0);
        Double_t value(0.0);
        
        if (arg < 0.0) {
                if (TMath::Abs(sigmaLVal_) > 1e-30) {
                        coef = -0.5/(sigmaLVal_*sigmaLVal_);
                }
        } else {
                if (TMath::Abs(sigmaRVal_) > 1e-30) {
                        coef = -0.5/(sigmaRVal_*sigmaRVal_);
                }
        }
        value = TMath::Exp(coef*arg*arg);

        // Calculate the norm
        Double_t xscaleL = LauConstants::root2*sigmaLVal_;
        Double_t xscaleR = LauConstants::root2*sigmaRVal_;
        
        Double_t integral(0.0);
        Double_t norm(0.0);
        Double_t result(0.0);


        if ( this->getMaxAbscissa() < meanVal_ ) {
                integral = sigmaLVal_ * ( TMath::Erf((this->getMaxAbscissa() - meanVal_)/xscaleL) - TMath::Erf((this->getMinAbscissa() - meanVal_)/xscaleL));
        } else if ( this->getMinAbscissa() > meanVal_ ) {
                integral = sigmaRVal_ * (TMath::Erf((this->getMaxAbscissa() - meanVal_)/xscaleR) - TMath::Erf((this->getMinAbscissa() - meanVal_)/xscaleR));
        } else {
                integral = sigmaRVal_*TMath::Erf((this->getMaxAbscissa() -meanVal_)/xscaleR) - sigmaLVal_*TMath::Erf((this->getMinAbscissa() - meanVal_)/xscaleL);
        }
        
        norm = LauConstants::rootPiBy2*integral;

        // the result
        result = value/norm;
        this->setUnNormPDFVal(result);
}

void LauDPDepBifurGaussPdf::scalePars_poly(Double_t dpPos)
{
        Int_t power = 1;
        for (vector<Double_t>::const_iterator iter = meanCoeffs_.begin(); iter != meanCoeffs_.end(); ++iter) {
                Double_t coeff = (*iter);
                meanVal_ += coeff * TMath::Power(dpPos,power);
                ++power;
        }

        power = 1;
        for (vector<Double_t>::const_iterator iter = sigmaLCoeffs_.begin(); iter != sigmaLCoeffs_.end(); ++iter) {
                Double_t coeff = (*iter);
                sigmaLVal_ += coeff * TMath::Power(dpPos,power);
                ++power;
        }

        power = 1;
        for (vector<Double_t>::const_iterator iter = sigmaRCoeffs_.begin(); iter != sigmaRCoeffs_.end(); ++iter) {
                Double_t coeff = (*iter);
                sigmaRVal_ += coeff * TMath::Power(dpPos,power);
                ++power;
        }
}

void LauDPDepBifurGaussPdf::scalePars_polyNegPower(Double_t dpPos)
{
        Int_t power = -1;
        for (vector<Double_t>::const_iterator iter = meanCoeffs_.begin(); iter != meanCoeffs_.end(); ++iter) {
                Double_t coeff = (*iter);
                meanVal_ += coeff * TMath::Power(dpPos,power);
                --power;
        }

        power = -1;
        for (vector<Double_t>::const_iterator iter = sigmaLCoeffs_.begin(); iter != sigmaLCoeffs_.end(); ++iter) {
                Double_t coeff = (*iter);
                sigmaLVal_ += coeff * TMath::Power(dpPos,power);
                --power;
        }

        power = -1;
        for (vector<Double_t>::const_iterator iter = sigmaRCoeffs_.begin(); iter != sigmaRCoeffs_.end(); ++iter) {
                Double_t coeff = (*iter);
                sigmaRVal_ += coeff * TMath::Power(dpPos,power);
                --power;
        }
}

void LauDPDepBifurGaussPdf::calcNorm() 
{
        this->setNorm(1.0);
}

void LauDPDepBifurGaussPdf::calcPDFHeight( const LauKinematics* kinematics )
{
        // Get the up to date parameter values
        meanVal_ = mean_->unblindValue();
        sigmaLVal_ = sigmaL_->unblindValue();
        sigmaRVal_ = sigmaR_->unblindValue();

        // Scale the gaussian parameters by the dpCentreDist (if appropriate)
        Double_t dpCentreDist = kinematics->distanceFromDPCentre();
        ScaleMethod scale = this->scaleMethod();
        if ( scale==poly ) {
                this->scalePars_poly(dpCentreDist);
        } else if ( scale==polyNegPower ) {
                this->scalePars_polyNegPower(dpCentreDist);
        } else {
                cerr<<"Scaling method unknown!  Methods available: <poly> and <polyNegPower>."<<endl;
                gSystem->Exit(EXIT_FAILURE);
        }
        
        // Calculate the PDF height for the Bifurcated Gaussian function.
        LauAbscissas maxPoint(3);
        maxPoint[0] = meanVal_;
        maxPoint[1] = kinematics->getm13Sq();
        maxPoint[2] = kinematics->getm23Sq();
        if ( meanVal_ > this->getMaxAbscissa() ) {
                maxPoint[0] = this->getMaxAbscissa();
        } else if ( meanVal_ < this->getMinAbscissa() ) {
                maxPoint[0] = this->getMinAbscissa();
        }       
        this->calcLikelihoodInfo(maxPoint);
        Double_t height = this->getUnNormLikelihood();

        // Multiply by a small factor to avoid problems from rounding errors
        height *= 1.001;

        this->setMaxHeight(height);
}