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

roofit/roofit/CMakeLists.txt

@@ -44,6 +44,7 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFit
     RooFunctor1DBinding.h
     RooFunctorBinding.h
     RooGamma.h
+    RooGaussExpTailDas.h
     RooGaussian.h
     RooGaussModel.h
     RooGExpModel.h
@@ -104,10 +105,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/RooGaussExpTailDas.cxx
     src/RooGaussian.cxx
     src/RooGaussModel.cxx
     src/RooGExpModel.cxx

roofit/roofit/inc/LinkDef1.h

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

roofit/roofit/inc/RooGaussExpTailDas.h

@@ -0,0 +1,48 @@
+#ifndef RooFit_RooFit_RooGaussExpTailDas_h
+#define RooFit_RooFit_RooGaussExpTailDas_h
+
+#include "RooAbsPdf.h"
+#include "RooRealProxy.h"
+
+class RooAbsReal;
+
+class RooGaussExpTailDas : public RooAbsPdf {
+public:
+  RooGaussExpTailDas() { };
+  RooGaussExpTailDas(const char *name, const char *title,
+                         RooAbsReal::Ref _x,
+                         RooAbsReal::Ref _x0,
+                         RooAbsReal::Ref _sigma,
+                         RooAbsReal::Ref _kL,
+                         RooAbsReal::Ref _kH
+                         );
+  RooGaussExpTailDas(const RooGaussExpTailDas& other, const char* name=nullptr);
+  TObject* clone(const char* newname) const override {
+    return new RooGaussExpTailDas(*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::Ref x() const { return *x_; }
+  RooAbsReal::Ref x0() const { return *x0_; }
+  RooAbsReal::Ref sigma() const { return *sigma_; }
+  RooAbsReal::Ref kL() const { return *kL_; }
+  RooAbsReal::Ref kH() const { return *kH_; }
+
+protected:
+  double evaluate() const override;
+
+private:
+  RooRealProxy x_;
+  RooRealProxy x0_;
+  RooRealProxy sigma_;
+  RooRealProxy kL_;
+  RooRealProxy kH_;
+
+private:
+
+  ClassDefOverride(RooGaussExpTailDas,1) // Gaussian with double-sided exponential tails PDF by DAS https://arxiv.org/abs/1603.08591v1
+};
+
+#endif

roofit/roofit/src/RooGaussExpTailDas.cxx

@@ -0,0 +1,110 @@
+/** \class RooGaussExpTailDas
+    \ingroup Roofit
+
+PDF implementing the Das 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
+\note To use the one-sided version, just set the opposite parameter k to a very large value
+*/
+
+
+#include "RooGaussExpTailDas.h"
+#include "RooAbsReal.h"
+#include <cmath>
+#include "Math/ProbFuncMathCore.h"
+
+ClassImp(RooGaussExpTailDas)
+
+
+//_____________________________________________________________________________
+RooGaussExpTailDas::RooGaussExpTailDas(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)
+{
+}
+
+
+//_____________________________________________________________________________
+RooGaussExpTailDas::RooGaussExpTailDas(const RooGaussExpTailDas& 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)
+{
+   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 RooGaussExpTailDas::evaluate() const
+{
+   Double_t t=(x_-x0_)/sigma_;
+
+   if (t<=-kL_)
+      return exp(kL_*kL_/2.+kL_*t);
+   else if (t>kH_)
+      return exp(kH_*kH_/2.-kH_*t);
+   else
+      return exp(-0.5*t*t);
+}
+
+
+//_____________________________________________________________________________
+Int_t RooGaussExpTailDas::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
+{
+   if( matchArgs(allVars,analVars,x_) )
+      return 1;
+
+   return 0;
+}
+
+
+//_____________________________________________________________________________
+Double_t RooGaussExpTailDas::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, TMath::Min(tmax, -kL_), kL_);
+   if (tmin <= kH_ && tmax > -kL_)
+      result += gaussianIntegral(TMath::Max(tmin, -kL_), TMath::Min(tmax, kH_));
+   if (tmax > kH_)
+      result += tailIntegral(TMath::Max(tmin, kH_), tmax, -kH_);
+
+   return sig*result;
+}