LauDPDepBifurGaussPdf.cc

Go to the documentation of this file.

/*
Copyright 2007 University of Warwick

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/

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