LauCrystalBallPdf.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

/*****************************************************************************
 * Class based on RooFit/RooCBShape.                                         *
 * Original copyright given below.                                           *
 *****************************************************************************
 * Authors:                                                                  *
 *   WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu       *
 *   DK, David Kirkby,    UC Irvine,         dkirkby@uci.edu                 *
 *                                                                           *
 * Copyright (c) 2000-2005, Regents of the University of California          *
 *                          and Stanford University. All rights reserved.    *
 *                                                                           *
 * Redistribution and use in source and binary forms,                        *
 * with or without modification, are permitted according to the terms        *
 * listed in LICENSE (http://roofit.sourceforge.net/license.txt)             *
 *****************************************************************************/

#include <iostream>
#include <vector>
using std::cout;
using std::cerr;
using std::endl;
using std::vector;

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

#include "LauConstants.hh"
#include "LauCrystalBallPdf.hh"

ClassImp(LauCrystalBallPdf)


LauCrystalBallPdf::LauCrystalBallPdf(const TString& theVarName, const vector<LauAbsRValue*>& params, Double_t minAbscissa, Double_t maxAbscissa) :
        LauAbsPdf(theVarName, params, minAbscissa, maxAbscissa),
        mean_(0),
        sigma_(0),
        alpha_(0),
        n_(0)
{
        // Constructor for the Crystal Ball PDF, which is a gaussian and a decaying tail
        // smoothly matched up. The tail goes as a 1/x^n
        //
        // The parameters in params are the mean and the sigma (half the width) of the gaussian,
        // the distance from the mean in which the gaussian and the tail are matched up (which
        // can be negative or positive), and the power "n" for the tail.
        // The last two arguments specify the range in which the PDF is defined, and the PDF
        // will be normalised w.r.t. these limits.

        mean_ = this->findParameter("mean");
        sigma_ = this->findParameter("sigma");
        alpha_ = this->findParameter("alpha");
        n_ = this->findParameter("order");

        if ((this->nParameters() != 4) || (mean_ == 0) || (sigma_ == 0) || (alpha_ == 0) || (n_ == 0)) {
                cerr<<"ERROR in LauCrystalBallPdf constructor: LauCrystalBallPdf requires 4 parameters: \"mean\", \"sigma\", \"alpha\" and \"order\"."<<endl;
                gSystem->Exit(EXIT_FAILURE);
        }

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

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

LauCrystalBallPdf::LauCrystalBallPdf(const LauCrystalBallPdf& other) : LauAbsPdf(other.varName(), other.getParameters(), other.getMinAbscissa(), other.getMaxAbscissa())
{
        // Copy constructor
        this->setRandomFun(other.getRandomFun());
        this->calcNorm();
}

void LauCrystalBallPdf::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];

        // Get the up to date parameter values
        Double_t mean = mean_->value();
        Double_t sigma = sigma_->value();
        Double_t alpha = alpha_->value();
        Double_t n = n_->value();

        Double_t result(0.0);
        Double_t t = (abscissa - mean)/sigma;
        if (alpha < 0.0) {
                t = -t;
        }

        Double_t absAlpha = TMath::Abs(alpha);

        if (t >= -absAlpha) {
                result = TMath::Exp(-0.5*t*t);
        } else {
                Double_t a = TMath::Power(n/absAlpha,n)*TMath::Exp(-0.5*absAlpha*absAlpha);
                Double_t b = n/absAlpha - absAlpha; 

                result = a/TMath::Power(b - t, n);
        }

        this->setUnNormPDFVal(result);
        
        // if the parameters are floating then we
        // need to recalculate the normalisation
        if (!this->cachePDF() && !this->withinNormCalc() && !this->withinGeneration()) {
                this->calcNorm();
        }

}

void LauCrystalBallPdf::calcNorm() 
{
        // Get the up to date parameter values
        Double_t mean = mean_->value();
        Double_t sigma = sigma_->value();
        Double_t alpha = alpha_->value();
        Double_t n = n_->value();

        Double_t result = 0.0;
        Bool_t useLog = kFALSE;

        if ( TMath::Abs(n-1.0) < 1.0e-05 ) {
                useLog = kTRUE;
        }

        Double_t sig = TMath::Abs(sigma);

        Double_t tmin = (this->getMinAbscissa() - mean)/sig;
        Double_t tmax = (this->getMaxAbscissa() - mean)/sig;
  
        if (alpha < 0) {
                Double_t tmp = tmin;
                tmin = -tmax;
                tmax = -tmp;
        }

        Double_t absAlpha = TMath::Abs(alpha);

        if ( tmin >= -absAlpha ) {
                result += sig*LauConstants::rootPiBy2*(   this->approxErf( tmax / LauConstants::root2 )
                                                        - approxErf( tmin / LauConstants::root2 ) );
        } else if ( tmax <= -absAlpha ) {
                Double_t a = TMath::Power(n/absAlpha, n)*TMath::Exp( -0.5*absAlpha*absAlpha);
                Double_t b = n/absAlpha - absAlpha;

                if ( useLog == kTRUE ) {
                        result += a*sig*( TMath::Log(b - tmin) - TMath::Log(b - tmax) );
                } else {
                        result += a*sig/(1.0 - n)*( 1.0/(TMath::Power( b - tmin, n - 1.0))
                                                - 1.0/(TMath::Power( b - tmax, n - 1.0)) );
                }
        } else {
                Double_t a = TMath::Power(n/absAlpha, n)*TMath::Exp( -0.5*absAlpha*absAlpha );
                Double_t b = n/absAlpha - absAlpha;

                Double_t term1 = 0.0;
                if ( useLog == kTRUE )
                        term1 = a*sig*( TMath::Log(b - tmin) - TMath::Log(n / absAlpha));
                else
                        term1 = a*sig/(1.0 - n)*( 1.0/(TMath::Power( b - tmin, n - 1.0))
                                                - 1.0/(TMath::Power( n/absAlpha, n - 1.0)) );
    
                Double_t term2 = sig*LauConstants::rootPiBy2*( this->approxErf( tmax / LauConstants::root2 )
                                     - this->approxErf( -absAlpha / LauConstants::root2 ) );

                result += term1 + term2;
        }
        this->setNorm(result);
}

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

        // Get the up to date parameter values
        Double_t mean = mean_->value();

        // The Crystall Ball function is a Gaussian with an exponentially decaying tail
        // Therefore, calculate the PDF height for the Gaussian function.
        LauAbscissas abscissa(1);
        abscissa[0] = mean;
        this->calcLikelihoodInfo(abscissa);

        Double_t height = this->getUnNormLikelihood();
        this->setMaxHeight(height);
}

Double_t LauCrystalBallPdf::approxErf(Double_t arg) const 
{
        static const Double_t erflim = 5.0;
        if ( arg > erflim ) {
                return 1.0;
        }
        if ( arg < -erflim ) {
                return -1.0;
        }

        return TMath::Erf(arg);
}