diff --git a/interface/AtlasPdfs.h b/interface/AtlasPdfs.h index 11964a8ba66..6aa9fe2ee36 100644 --- a/interface/AtlasPdfs.h +++ b/interface/AtlasPdfs.h @@ -388,3 +388,100 @@ class RooStarMomentMorph : public RooAbsPdf { #endif +// --- Importing file RooTwoSidedCBShape +#ifndef ROOT_Hfitter_RooTwoSidedCBShape +#define ROOT_Hfitter_RooTwoSidedCBShape + +#include "RooAbsPdf.h" +#include "RooRealProxy.h" +#include "RooAbsReal.h" +#include "RooArgSet.h" + +class RooRealVar; + + +class RooTwoSidedCBShape : public RooAbsPdf { +public: + RooTwoSidedCBShape() {} + RooTwoSidedCBShape(const char *name, const char *title, RooAbsReal& _m, + RooAbsReal& _m0, RooAbsReal& _sigma, + RooAbsReal& _alphaLo, RooAbsReal& _nLo, + RooAbsReal& _alphaHi, RooAbsReal& _nHi); + + RooTwoSidedCBShape(const RooTwoSidedCBShape& other, const char* name = 0); + virtual TObject* clone(const char* newname) const { return new RooTwoSidedCBShape(*this,newname); } + + inline virtual ~RooTwoSidedCBShape() { } + + Int_t getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName=0) const; + Double_t analyticalIntegral(Int_t code, const char* rangeName=0) const; + +protected: + + Double_t gaussianIntegral(Double_t tmin, Double_t tmax) const; + Double_t powerLawIntegral(Double_t tmin, Double_t tmax, Double_t alpha, Double_t n) const; + Double_t evaluate() const; + + RooRealProxy m; + RooRealProxy m0; + RooRealProxy sigma; + RooRealProxy alphaLo; + RooRealProxy nLo; + RooRealProxy alphaHi; + RooRealProxy nHi; + + +private: + + ClassDef(RooTwoSidedCBShape,1) +}; + + +#endif + +// --- Importing file HggBernstein +/***************************************************************************** + * Authors: * + * Nicolas Berger (nicolas.berger@cern.ch) + * * + * * + * 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) * + *****************************************************************************/ +#ifndef ROOT_HggBernstein +#define ROOT_HggBernstein + +#include "RooAbsPdf.h" +#include "RooRealProxy.h" +#include "RooListProxy.h" + +class RooRealVar; +class RooArgList ; + +class HggBernstein : public RooAbsPdf { + public: + + HggBernstein() ; + HggBernstein(const char *name, const char *title, + RooAbsReal& _x, const RooArgList& _coefList, double xMin = 0, double xMax = -1) ; + + HggBernstein(const HggBernstein& other, const char* name = 0); + virtual TObject* clone(const char* newname) const { return new HggBernstein(*this, newname); } + inline virtual ~HggBernstein() { } + + Int_t getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName=0) const ; + Double_t analyticalIntegral(Int_t code, const char* rangeName=0) const ; + + private: + + RooRealProxy _x; + RooListProxy _coefList ; + double _xMin, _xMax; + + Double_t evaluate() const; + + ClassDef(HggBernstein,1) // Bernstein polynomial PDF + }; +#endif + diff --git a/src/AtlasPdfs.cxx b/src/AtlasPdfs.cxx index 34e758a3ad2..8967190c326 100644 --- a/src/AtlasPdfs.cxx +++ b/src/AtlasPdfs.cxx @@ -1924,3 +1924,226 @@ Bool_t RooStarMomentMorph::setBinIntegrator(RooArgSet& allVars) } return false; } + +// --- Importing file RooTwoSidedCBShape +#include "RooFit.h" + +#include "Riostream.h" +#include "Riostream.h" +#include + +#include "RooAbsReal.h" +#include "RooRealVar.h" +#include "RooMath.h" +#include "TMath.h" +#include "Math/ProbFuncMathCore.h" + +ClassImp(RooTwoSidedCBShape) + +//_____________________________________________________________________________ +RooTwoSidedCBShape::RooTwoSidedCBShape(const char *name, const char *title, + RooAbsReal& _m, RooAbsReal& _m0, RooAbsReal& _sigma, + RooAbsReal& _alphaLo, RooAbsReal& _nLo, + RooAbsReal& _alphaHi, RooAbsReal& _nHi) : + RooAbsPdf(name, title), + m("m", "Dependent", this, _m), + m0("m0", "M0", this, _m0), + sigma("sigma", "Sigma", this, _sigma), + alphaLo("alphaLo", "Low-side Alpha", this, _alphaLo), + nLo("nLo", "Low-side Order", this, _nLo), + alphaHi("alphaHi", "High-side Alpha", this, _alphaHi), + nHi("nHi", "Hig-side Order", this, _nHi) +{ +} + + +//_____________________________________________________________________________ +RooTwoSidedCBShape::RooTwoSidedCBShape(const RooTwoSidedCBShape& other, const char* name) : + RooAbsPdf(other, name), m("m", this, other.m), m0("m0", this, other.m0), + sigma("sigma", this, other.sigma), + alphaLo("alphaLo", this, other.alphaLo), nLo("nLo", this, other.nLo), + alphaHi("alphaHi", this, other.alphaHi), nHi("nHi", this, other.nHi) +{ +} + + +//_____________________________________________________________________________ +Double_t RooTwoSidedCBShape::evaluate() const { + + Double_t t = (m-m0)/sigma; + + if (t < -alphaLo) { + Double_t a = exp(-0.5*alphaLo*alphaLo); + Double_t b = nLo/alphaLo - alphaLo; + return a/TMath::Power(alphaLo/nLo*(b - t), nLo); + } + else if (t > alphaHi) { + Double_t a = exp(-0.5*alphaHi*alphaHi); + Double_t b = nHi/alphaHi - alphaHi; + return a/TMath::Power(alphaHi/nHi*(b + t), nHi); + } + return exp(-0.5*t*t); +} + + +//_____________________________________________________________________________ +Int_t RooTwoSidedCBShape::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const +{ + if (matchArgs(allVars,analVars,m)) return 1; + return 0; +} + + +//_____________________________________________________________________________ +Double_t RooTwoSidedCBShape::analyticalIntegral(Int_t code, const char* rangeName) const +{ + assert(code==1); + Double_t result = 0; + + Double_t sig = fabs((Double_t)sigma); + Double_t tmin = (m.min(rangeName)-m0)/sig; + Double_t tmax = (m.max(rangeName)-m0)/sig; + + if (tmin < -alphaLo) + result += powerLawIntegral(tmin, TMath::Min(tmax, -alphaLo), alphaLo, nLo); + if (tmin < alphaHi && tmax > -alphaLo) + result += gaussianIntegral(TMath::Max(tmin, -alphaLo), TMath::Min(tmax, alphaHi)); + if (tmax > alphaHi) + result += powerLawIntegral(-tmax, TMath::Min(-tmin, -alphaHi), alphaHi, nHi); + + return sig*result; +} + + +//_____________________________________________________________________________ +Double_t RooTwoSidedCBShape::gaussianIntegral(Double_t tmin, Double_t tmax) const +{ + return sqrt(TMath::TwoPi())*(ROOT::Math::gaussian_cdf(tmax) - ROOT::Math::gaussian_cdf(tmin)); +} + + +//_____________________________________________________________________________ +Double_t RooTwoSidedCBShape::powerLawIntegral(Double_t tmin, Double_t tmax, Double_t alpha, Double_t n) const +{ + // Implement protection for n = 1 from RooCBShape.cxx + bool useLog = false; + if( fabs(n-1.0) < 1.0e-05 ) useLog = true; + + Double_t a = exp(-0.5*alpha*alpha); + Double_t b = n/alpha - alpha; + + if( useLog ) return a * TMath::Power(n/alpha, n) * ( log( b-tmin ) - log( b-tmax ) ); + return a/(1 - n)*( (b - tmin)/(TMath::Power(alpha/n*(b - tmin), n)) - (b - tmax)/(TMath::Power(alpha/n*(b - tmax), n)) ); +} + +// --- Importing file HggBernstein + +#include "RooFit.h" +#include "Riostream.h" +#include "Riostream.h" +#include +#include "TMath.h" +#include "RooAbsReal.h" +#include "RooRealVar.h" +#include "RooArgList.h" + +#include +using std::cout; +using std::endl; + +ClassImp(HggBernstein) +//_____________________________________________________________________________ +HggBernstein::HggBernstein() :_xMin(0), _xMax(1) +{ +} + + +//_____________________________________________________________________________ +HggBernstein::HggBernstein(const char* name, const char* title, + RooAbsReal& x, const RooArgList& coefList, double xMin, double xMax): + RooAbsPdf(name, title), + _x("x", "Dependent", this, x), + _coefList("coefficients","List of coefficients",this), + _xMin(xMin < xMax ? xMin : _x.min()), + _xMax(xMin < xMax ? xMax : _x.max()) +{ + //cout << _xMin << " " << _xMax << " " << xMin << " " << xMax << " " << _x.min() << " " << _x.max() <Next())) { + if (!dynamic_cast(coef)) { + cout << "HggBernstein::ctor(" << GetName() << ") ERROR: coefficient " << coef->GetName() + << " is not of type RooAbsReal" << endl ; + assert(0) ; + } + _coefList.add(*coef) ; + } + delete coefIter ; +} + + + +//_____________________________________________________________________________ +HggBernstein::HggBernstein(const HggBernstein& other, const char* name) : + RooAbsPdf(other, name), + _x("x", this, other._x), + _coefList("coefList",this,other._coefList), + _xMin(other._xMin), _xMax(other._xMax) +{ +} + + +//_____________________________________________________________________________ +Double_t HggBernstein::evaluate() const +{ + Int_t degree= _coefList.getSize()-1; // n+1 polys of degree n + + Double_t temp=0, tempx = (_x-_xMin)/(_xMax-_xMin); // rescale to [0,1] + + RooFIter iter = _coefList.fwdIterator() ; + for (int i=0; i<=degree; ++i){ + temp += ((RooAbsReal*)iter.next())->getVal() * + TMath::Binomial(degree, i) * pow(tempx,i) * pow(1-tempx,degree-i); + } + return temp; + +} + + +//_____________________________________________________________________________ +Int_t HggBernstein::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const +{ + // No analytical calculation available (yet) of integrals over subranges + if (rangeName && strlen(rangeName)) { + return 0 ; + } + + if (matchArgs(allVars, analVars, _x)) return 1; + return 0; +} + + +//_____________________________________________________________________________ +Double_t HggBernstein::analyticalIntegral(Int_t code, const char* rangeName) const +{ + assert(code==1) ; + Double_t xmin = _x.min(rangeName); Double_t xmax = _x.max(rangeName); + Double_t umin = (xmin - _xMin)/(_xMax - _xMin); + Double_t umax = (xmax - _xMin)/(_xMax - _xMin); + Int_t degree= _coefList.getSize()-1; // n+1 polys of degree n + Double_t norm(0) ; + + // for each of the i Bernstein basis polynomials + // represent it in the 'power basis' (the naive polynomial basis) + // where the integral is straight forward. + for (int j = 0; j <= degree; ++j){ // power basis + Double_t sum = 0; + for (int i = 0; i <= j; ++i) sum += ((j-i) % 2 ? -1 : 1)*TMath::Binomial(j,i)*((RooAbsReal*)_coefList.at(i))->getVal(); + sum *= TMath::Binomial(degree, j); + sum *= (TMath::Power(umax, j+1) - TMath::Power(umin, j+1))/(j+1); + norm += sum; + } + norm *= (_xMax - _xMin); + return norm; +} diff --git a/src/classes_def.xml b/src/classes_def.xml index ae303590427..92cce488fdd 100644 --- a/src/classes_def.xml +++ b/src/classes_def.xml @@ -172,6 +172,8 @@ + +