[RF] New PDF: gaussian with double sided exponential tails

README/ReleaseNotes/v638/index.md

@@ -47,6 +47,11 @@ RooAbsData *dataForChan = found != splits.end() ? found->get() : nullptr;
 // ... do something with channelData ...
 ```
 
+### RooCrystalBall alternative
+
+- A [simpler alternative](http://arxiv.org/abs/1603.08591v1) to RooCrystalBall using a Gaussian with exponential tails has been implemented: `class RooGaussExpTails`.
+
+
 ## RDataFrame
 - Memory savings in RDataFrame: When many Histo3D are filled in RDataFrame, the memory consumption in multi-threaded runs can be prohibitively large, because
   RDF uses one copy of each histogram per thread. Now, RDataFrame can reduce the number of clones using `ROOT::RDF::Experimental::ThreadsPerTH3()`. Setting this

roofit/roofit/CMakeLists.txt

@@ -45,6 +45,7 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFit
     RooFunctor1DBinding.h
     RooFunctorBinding.h
     RooGamma.h
+    RooGaussExpTails.h
     RooGaussian.h
     RooGaussModel.h
     RooGExpModel.h
@@ -106,10 +107,11 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFit
     src/RooDstD0BG.cxx
     src/RooExponential.cxx
     src/RooLegacyExpPoly.cxx
-        src/RooPowerSum.cxx
+    src/RooPowerSum.cxx
     src/RooFunctor1DBinding.cxx
     src/RooFunctorBinding.cxx
     src/RooGamma.cxx
+    src/RooGaussExpTails.cxx
     src/RooGaussian.cxx
     src/RooGaussModel.cxx
     src/RooGExpModel.cxx

roofit/roofit/inc/LinkDef1.h

@@ -23,6 +23,7 @@
 #pragma link C++ class RooExponential+ ;
 #pragma link C++ class RooLegacyExpPoly+ ;
 #pragma link C++ class RooPowerSum+ ;
+#pragma link C++ class RooGaussExpTails+ ;
 #pragma link C++ class RooGaussian+ ;
 #pragma link C++ class RooLognormal+ ;
 #pragma link C++ class RooGamma+ ;

roofit/roofit/inc/RooGaussExpTails.h

@@ -0,0 +1,40 @@
+#ifndef RooFit_RooFit_RooGaussExpTails_h
+#define RooFit_RooFit_RooGaussExpTails_h
+
+#include "RooAbsPdf.h"
+#include "RooRealProxy.h"
+
+class RooAbsReal;
+
+class RooGaussExpTails : public RooAbsPdf {
+public:
+   RooGaussExpTails() {}
+   RooGaussExpTails(const char *name, const char *title, RooAbsReal::Ref x, RooAbsReal::Ref x0, RooAbsReal::Ref sigma,
+                    RooAbsReal::Ref kL, RooAbsReal::Ref kH);
+   RooGaussExpTails(const RooGaussExpTails &other, const char *name = nullptr);
+   TObject *clone(const char *newname) const override { return new RooGaussExpTails(*this, newname); }
+
+   Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName = nullptr) const override;
+   double analyticalIntegral(Int_t code, const char *rangeName = nullptr) const override;
+
+   RooAbsReal const &x() const { return *_x; }
+   RooAbsReal const &x0() const { return *_x0; }
+   RooAbsReal const &sigma() const { return *_sigma; }
+   RooAbsReal const &kL() const { return *_kL; }
+   RooAbsReal const &kH() const { return *_kH; }
+
+protected:
+   double evaluate() const override;
+
+private:
+   RooRealProxy _x;
+   RooRealProxy _x0;
+   RooRealProxy _sigma;
+   RooRealProxy _kL;
+   RooRealProxy _kH;
+
+private:
+   ClassDefOverride(RooGaussExpTails, 1) // Gaussian with double-sided exponential tails PDF, see https://arxiv.org/abs/1603.08591v1
+};
+
+#endif

roofit/roofit/src/RooGaussExpTails.cxx

@@ -0,0 +1,98 @@
+/** \class RooGaussExpTails
+    \ingroup Roofit
+
+PDF implementing a Gaussian core + double-sided exponential tail distribution
+\author Souvik Das (8/1/2013) Initial implementation and
+Giovanni Marchiori (30/3/2016) Implemented analytic integral
+\see http://arxiv.org/pdf/1603.08591v1.pdf, https://github.com/mpuccio/RooCustomPdfs/blob/master/RooGausDExp.cxx, https://doi.org/10.1142/S0217751X14300440
+\note To use the one-sided version, just set the opposite parameter k to a very large value
+*/
+
+
+#include "RooGaussExpTails.h"
+#include "RooAbsReal.h"
+#include <cmath>
+#include "Math/ProbFuncMathCore.h"
+
+//_____________________________________________________________________________
+RooGaussExpTails::RooGaussExpTails(const char *name, const char *title, RooAbsReal::Ref x, RooAbsReal::Ref x0,
+                                   RooAbsReal::Ref sigma, RooAbsReal::Ref kL, RooAbsReal::Ref kH)
+   : RooAbsPdf(name, title),
+     _x("x", "x", this, x),
+     _x0("x0", "x0", this, x0),
+     _sigma("sigma", "sigma", this, sigma),
+     _kL("kL", "kL", this, kL),
+     _kH("kH", "kH", this, kH)
+{
+}
+
+//_____________________________________________________________________________
+RooGaussExpTails::RooGaussExpTails(const RooGaussExpTails &other, const char* name)
+   : RooAbsPdf(other, name),
+     _x("x", this, other._x),
+     _x0("x0", this, other._x0),
+     _sigma("sigma", this, other._sigma),
+     _kL("kL", this, other._kL),
+     _kH("kH", this, other._kH)
+{
+}
+
+////////////////////////////////////////////////////////////////////////////////
+
+namespace {
+
+inline double gaussianIntegral(double tmin, double tmax)
+{
+   constexpr double m_sqrt_2_pi = 2.50662827463; // std::sqrt(TMath::TwoPi())
+   return m_sqrt_2_pi * (ROOT::Math::gaussian_cdf(tmax) - ROOT::Math::gaussian_cdf(tmin));
+}
+
+inline double tailIntegral(double tmin, double tmax, double k)
+{
+   double a = std::exp(0.5 * k * k) / k;
+   return a * (std::exp(k * tmax) - std::exp(k * tmin));
+}
+
+} // namespace
+
+//_____________________________________________________________________________
+Double_t RooGaussExpTails::evaluate() const
+{
+   Double_t t = (_x - _x0) / _sigma;
+
+   if (t <= -_kL)
+      return std::exp(0.5 * _kL * _kL + _kL * t);
+   else if (t > _kH)
+      return std::exp(0.5 * _kH * _kH - _kH * t);
+   else
+      return std::exp(-0.5 * t * t);
+}
+
+//_____________________________________________________________________________
+Int_t RooGaussExpTails::getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char * /*rangeName*/) const
+{
+   if (matchArgs(allVars, analVars, _x))
+      return 1;
+
+   return 0;
+}
+
+//_____________________________________________________________________________
+Double_t RooGaussExpTails::analyticalIntegral(Int_t code, const char *rangeName) const
+{
+   R__ASSERT(code == 1);
+   double result = 0;
+
+   double sig = std::abs((Double_t)_sigma);
+   double tmin = (_x.min(rangeName) - _x0) / sig;
+   double tmax = (_x.max(rangeName) - _x0) / sig;
+
+   if (tmin <= -_kL)
+      result += tailIntegral(tmin, std::min(tmax, -_kL), _kL);
+   if (tmin <= _kH && tmax > -_kL)
+      result += gaussianIntegral(std::max(tmin, -_kL), std::min(tmax, +_kH));
+   if (tmax > _kH)
+      result += tailIntegral(std::max(tmin, +_kH), tmax, -_kH);
+
+   return sig * result;
+}