From 624ec497835ba644869c2c3e0c45d4ce1318ae2f Mon Sep 17 00:00:00 2001 From: Vanya Belyaev Date: Thu, 21 Sep 2023 12:27:55 +0200 Subject: [PATCH] 1. Add `Ostap::MoreRooFit::Rational` and `Ostap::MoreRooFit::RatioalBernstein` 1. Add `RationaFun` FUN --- ReleaseNotes/release_notes.md | 3 + ostap/fitting/roofuncs.py | 66 ++++++ ostap/fitting/rooreduce.py | 8 +- ostap/fitting/tests/test_fitting_transform.py | 20 +- source/include/Ostap/MoreRooFit.h | 22 ++ source/include/Ostap/MoreVars.h | 173 ++++++++++++++- source/include/Ostap/Rational.h | 2 +- source/src/MoreRooFit.cpp | 33 ++- source/src/MoreVars.cpp | 198 +++++++++++++++++- 9 files changed, 504 insertions(+), 21 deletions(-) diff --git a/ReleaseNotes/release_notes.md b/ReleaseNotes/release_notes.md index 8bd19a20..5d740306 100644 --- a/ReleaseNotes/release_notes.md +++ b/ReleaseNotes/release_notes.md @@ -1,5 +1,8 @@ ## New features: + 1. Add `Ostap::MoreRooFit::Rational` and `Ostap::MoreRooFit::RatioalBernstein` + 1. Add `RationaFun` FUN + ## Backward incompatible: ## Bug fixes: diff --git a/ostap/fitting/roofuncs.py b/ostap/fitting/roofuncs.py index 14585d1b..7aae29ac 100644 --- a/ostap/fitting/roofuncs.py +++ b/ostap/fitting/roofuncs.py @@ -21,6 +21,7 @@ 'RooPoly' , ## simple wrapper for RooPolyVar (RooAbsReal) 'ScaleAndShift' , ## scale and shift (RooAbsReal) 'BSplineFun' , ## BSpline (RooAbsReal) + 'RationalFun' , ## Ratioal function (RooAbsReal) 'Shape1D_fun' , ## arbitrary fixed shape (RooAbsReal) 'Histo1D_fun' , ## fixed shap form historgam (RooAbsReal) 'Histo1DErr_fun' , ## fixed shap from historgam errors (RooAbsReal) @@ -538,6 +539,71 @@ def power ( self ) : """'power' : degree for BSpline object""" return self.fun.degree() + + +# ============================================================================= +## @class RationalFun +# A simple pole-free rational function at interval \f$ x_{min} \le x \le x_{max}\f$ +# \f[ F(x) = \frac{p(x)}{q(x)} \f] +# Actually internally it uses +# the Floater-Hormann rational barycentric interpolant +# and parameters are the function values at Chebyshev's nodes +# +# @see Ostap::MoreRooFit::Rational +# @see Ostap::Math::Rational +# @see Ostap::Math::FloaterHormann +# @author Vanya BELYAEV Ivan.Belyaev@cern.ch +# @date 2023-09-21 +class RationalFun(FUN1,ParamsPoly) : + r"""A simple pole-free rational function at interval \f$ x_{min} \le x \le x_{max}\f$ + >>> p1 = BernsteinPoly ( 'P1' , xvar = (0,1) , power = 3 ) + >>> p2 = BernsteinPoly ( 'P2' , xvar = (0,1) , pars = ... ) + - see Ostap.MoreRooFit.Bernstein + - see Ostap.Math.Bernstein + """ + def __init__ ( self , name , xvar , + n = 3 , + d = 1 , pars = None ) : + + ## initialize the base class + FUN1 .__init__ ( self , name , xvar = xvar ) + ParamsPoly.__init__ ( self , + npars = n + 1 , + pars = pars ) + + assert isinstance ( n , int ) and 1 <= n , 'Invalid parameter n' + assert isinstance ( d , int ) and 0 <= d <= n , 'Invalid parameter d' + + xmin , xmax = self.xminmax () + + ## create the function + self.fun = Ostap.MoreRooFit.Rational ( + self.roo_name ( 'rfun_' ) , + 'Rational %s' % self.name , + self.xvar , + self.pars_lst , + d , + xmin , + xmax ) + + ## self.tricks = True + self.config = { + 'name' : self.name , + 'xvar' : self.xvar , + 'n' : self.fun.n () , + 'd' : self.fun.d () , + 'pars' : self.pars , + } + + @property + def n ( self ) : + """'n' : n-parameter of rational function""" + return self.fun.n () + @property + def d ( self ) : + """'d' : d-parameter of rational function""" + return self.fun.d () + # ============================================================================= ## Generic 1D-shape from C++ callable diff --git a/ostap/fitting/rooreduce.py b/ostap/fitting/rooreduce.py index 55d13f45..8fa59b88 100644 --- a/ostap/fitting/rooreduce.py +++ b/ostap/fitting/rooreduce.py @@ -2668,7 +2668,8 @@ def _rgauss3d_reduce_ ( pdf ): # ============================================================================= ## reduce Rational def _rational_reduce_ ( pdf ): - """Reduce Ostap::Models::Rational""" + """Reduce Ostap::Models::Rational & Ostap::MoreRooFit::Rational + """ return root_store_factory , ( type ( pdf ) , pdf.name , pdf.title , @@ -2678,7 +2679,10 @@ def _rational_reduce_ ( pdf ): pdf.xmin () , pdf.xmax () ) -Ostap.Models.Rational.__reduce__ = _rational_reduce_ +Ostap.Models.Rational .__reduce__ = _rational_reduce_ +Ostap.MoreRooFit.Rational .__reduce__ = _rational_reduce_ +Ostap.MoreRooFit.RationalBernstein .__reduce__ = _rational_reduce_ + # ============================================================================= diff --git a/ostap/fitting/tests/test_fitting_transform.py b/ostap/fitting/tests/test_fitting_transform.py index 96dfd7a3..ec642330 100755 --- a/ostap/fitting/tests/test_fitting_transform.py +++ b/ostap/fitting/tests/test_fitting_transform.py @@ -14,7 +14,7 @@ __author__ = "Ostap developers" __all__ = () ## nothing to import # ============================================================================= -import ostap.fitting.roofit +from ostap.core.meta_info import root_info from ostap.core.core import cpp, VE, dsID import ostap.fitting.models as Models from ostap.fitting.transform_pdf import TrPDF @@ -23,6 +23,7 @@ from ostap.utils.timing import timing from ostap.plotting.canvas import use_canvas from ostap.utils.utils import wait +import ostap.fitting.roofit import ROOT, random # ============================================================================= # logging @@ -40,11 +41,14 @@ def test_transform () : logger = getLogger ( 'test_transform' ) - if not 62000 <= ROOT.gROOT.GetVersionInt() : - logger.info ("Not for this version of ROOT") + if not (6,20) <= root_info : + logger.info ("Not for this version of ROOT!") return + + lmin = 0 + lmax = 5 - x = ROOT.RooRealVar('x','',1,100000) + x = ROOT.RooRealVar('x','', 10.**lmin , 10.**lmax ) mean = 1000 sigma = 1000 @@ -58,15 +62,15 @@ def test_transform () : dataset.add ( varset ) dataset.add_var ( 'lx' , 'log10(x)' ) - dataset.lx.setMin ( 0 ) - dataset.lx.setMax ( 5 ) + dataset.lx.setMin ( lmin ) + dataset.lx.setMax ( lmax ) ## original PDF gauss = Models.Gauss_pdf ( 'G' , xvar = x , mean = ( mean , mean / 2 , mean * 2 ) , sigma = ( sigma , sigma / 2 , sigma * 2 ) ) - with use_canvas ( 'test_transform' ) , wait ( 1 ) : + with use_canvas ( 'test_transform/x' ) , wait ( 1 ) : r1 , f1 = gauss.fitTo ( dataset , draw = True , silent = True ) logger.info ( 'Fit x:\n%s' % r1.table() ) @@ -77,7 +81,7 @@ def test_transform () : ## transformed PDF tgauss = TrPDF ( pdf = gauss , new_var = NX ) - with use_canvas ( 'test_transform' ) , wait ( 1 ): + with use_canvas ( 'test_transform/log10' ) , wait ( 1 ): r2 , f2 = tgauss.fitTo ( dataset , draw = True , silent = True ) logger.info ( 'Fit log10(x):\n%s' % r2.table() ) diff --git a/source/include/Ostap/MoreRooFit.h b/source/include/Ostap/MoreRooFit.h index db429e93..9b736afd 100644 --- a/source/include/Ostap/MoreRooFit.h +++ b/source/include/Ostap/MoreRooFit.h @@ -2729,6 +2729,28 @@ namespace Ostap /// clone AddDeps* clone ( const char* newname ) const override ; // ====================================================================== + public: // delegation + // ====================================================================== + Double_t analyticalIntegral + ( Int_t code , + const char* range = nullptr ) const override ; + // + Double_t analyticalIntegralWN + ( Int_t code , + const RooArgSet* normset , + const char* range = nullptr ) const override ; + // + Int_t getAnalyticalIntegral + ( RooArgSet& allVars , + RooArgSet& analVars , + const char* range = nullptr ) const override ; + // + Int_t getAnalyticalIntegralWN + ( RooArgSet& allVars , + RooArgSet& analVars , + const RooArgSet* normset , + const char* range = nullptr ) const override ; + // ====================================================================== public: // ====================================================================== /// the actual evaluation of the result diff --git a/source/include/Ostap/MoreVars.h b/source/include/Ostap/MoreVars.h index 7852d414..7a1bf101 100644 --- a/source/include/Ostap/MoreVars.h +++ b/source/include/Ostap/MoreVars.h @@ -20,6 +20,7 @@ #include "Ostap/Bernstein.h" #include "Ostap/Bernstein1D.h" #include "Ostap/BSpline.h" +#include "Ostap/Rational.h" #include "Ostap/HistoInterpolators.h" // ============================================================================ /// forward declarations @@ -494,7 +495,177 @@ namespace Ostap // ====================================================================== } ; // The end of class Ostap::MoreRooFit::BSpline // ======================================================================== - /** @class Shape + /** @class Rational + * A simple pole-free rational function at interval \f$ x_{min} \le x \le x_{max}\f$ + * \f[ F(x) = \frac{p(x)}{q(x)} \f] + * Actually internally it uses + * the Floater-Hormann rational barycentric interpolant + * and parameters are the function valeus at Chebyshev's nodes + * + * @see Ostap::Math::Rational + * @see Ostap::Math::FloaterHormann + * @author Vanya BELYAEV Ivan.Belyaev@cern.ch + * @date 2023-09-21 + */ + class Rational final : public RooAbsReal + { + public: + // ====================================================================== + ClassDefOverride(Ostap::MoreRooFit::Rational, 1) ; + // ====================================================================== + public: + // ====================================================================== + Rational + ( const std::string& name , + const std::string& title , + RooAbsReal& xvar , + const RooArgList& pars , + const unsigned short d , + const double xmin , + const double xmax ) ; + /// copy constructor + Rational ( const Rational& right , const char* name = nullptr ) ; + /// default + Rational () ; + /// clone method + Rational* clone ( const char* name ) const override ; + /// virtual destructor + virtual ~Rational() ; + // ====================================================================== + public: + // ====================================================================== + Int_t getAnalyticalIntegral + ( RooArgSet& allVars , + RooArgSet& analVars , + const char* rangeName = nullptr ) const override ; + Double_t analyticalIntegral + ( Int_t code , + const char* rangeName = nullptr ) const override ; + // ====================================================================== + public: + // ====================================================================== + /// get the variable/observable + const RooAbsReal& xvar () const { return m_xvar.arg() ; } + /// get parameters + const RooArgList& pars () const { return m_pars ; } + /// get n + unsigned short n () const { return m_rational.n () ; } + /// get d + unsigned short d () const { return m_rational.d () ; } + /// xmin for Rational + double xmin () const { return m_rational.xmin () ; } + /// xmax for Rational + double xmax () const { return m_rational.xmax () ; } + // ====================================================================== + public: + // ====================================================================== + void setPars () const ; + // ====================================================================== + public: + // ====================================================================== + const Ostap::Math::Rational& function () const { return m_rational ; } + const Ostap::Math::Rational& rational () const { return m_rational ; } + // ====================================================================== + public: + // ====================================================================== + Double_t evaluate () const override ; + // ====================================================================== + private: + // ====================================================================== + RooRealProxy m_xvar {} ; + RooListProxy m_pars {} ; + mutable Ostap::Math::Rational m_rational { 3 , 1 } ; + // ====================================================================== + } ; + // ======================================================================== + /** @class RationalBernstein + * Rational fnuction as ratio of Bernstein polyhomial and + * positive Bernstein polynomial + * \f[ R ( x ) = \frac{B(x)}{P(x)\f] + * @see Ostap::Math::RationalBernstein + * @see Ostap::Math::Bernstein + * @see Ostap::Math::Positive + * @author Vanya BELYAEV Ivan.Belyaev@cern.ch + * @date 2023-09-21 + */ + class RationalBernstein final : public RooAbsReal + { + public: + // ====================================================================== + ClassDefOverride(Ostap::MoreRooFit::RationalBernstein, 1) ; + // ====================================================================== + public: + // ====================================================================== + RationalBernstein + ( const std::string& name , + const std::string& title , + RooAbsReal& xvar , + const RooArgList& p , // numerator + const RooArgList& q , // denominator + const double xmin , + const double xmax ) ; + RationalBernstein + ( const std::string& name , + const std::string& title , + RooAbsReal& xvar , + const RooArgList& pars , // all pars + const unsigned short p , // degree of numerator + const double xmin , + const double xmax ) ; + /// copy constructor + RationalBernstein ( const RationalBernstein& right , const char* name = nullptr ) ; + /// default + RationalBernstein () ; + /// clone method + RationalBernstein* clone ( const char* name ) const override ; + /// virtual destructor + virtual ~RationalBernstein() ; + // ====================================================================== + public: + // ====================================================================== + Int_t getAnalyticalIntegral + ( RooArgSet& allVars , + RooArgSet& analVars , + const char* rangeName = nullptr ) const override ; + Double_t analyticalIntegral + ( Int_t code , + const char* rangeName = nullptr ) const override ; + // ====================================================================== + public: + // ====================================================================== + /// get the variable/observable + const RooAbsReal& xvar () const { return m_xvar.arg() ; } + /// get parameters + const RooArgList& pars () const { return m_pars ; } + /// get p + unsigned short p () const { return m_rational.pdegree () ; } + /// xmin for Rational + double xmin () const { return m_rational.xmin () ; } + /// xmax for Rational + double xmax () const { return m_rational.xmax () ; } + // ====================================================================== + public: + // ====================================================================== + void setPars () const ; + // ====================================================================== + public: + // ====================================================================== + const Ostap::Math::RationalBernstein& function () const { return m_rational ; } + const Ostap::Math::RationalBernstein& rational () const { return m_rational ; } + // ====================================================================== + public: + // ====================================================================== + Double_t evaluate () const override ; + // ====================================================================== + private: + // ====================================================================== + RooRealProxy m_xvar {} ; + RooListProxy m_pars {} ; + mutable Ostap::Math::RationalBernstein m_rational { 3 , 1 } ; + // ====================================================================== + } ; + // ======================================================================== + /** @class Shape1D * The generic "fixed shape" function * @author Vanya BELYAEV Ivan.Belyaev@itep.ru * @date 2023-01-30 diff --git a/source/include/Ostap/Rational.h b/source/include/Ostap/Rational.h index e4a5cfd8..5da6ed9e 100644 --- a/source/include/Ostap/Rational.h +++ b/source/include/Ostap/Rational.h @@ -26,7 +26,7 @@ namespace Ostap * the Floater-Hormann rational barycentric interpolant * and parameters are the function valeus at Chebyshev's nodes * - # @see Ostap::Math::FloaterHormann + * @see Ostap::Math::FloaterHormann * @author Vanya BELYAEV Ivan.Belyaev@cern.ch * @date 2023-09-14 */ diff --git a/source/src/MoreRooFit.cpp b/source/src/MoreRooFit.cpp index 25b3259a..e154065b 100644 --- a/source/src/MoreRooFit.cpp +++ b/source/src/MoreRooFit.cpp @@ -1086,19 +1086,19 @@ Double_t Ostap::MoreRooFit::Id::analyticalIntegral ( Int_t code , const char* range ) const { return m_x.arg().analyticalIntegral ( code , range ) ; } -// +// ============================================================================ Double_t Ostap::MoreRooFit::Id::analyticalIntegralWN ( Int_t code , const RooArgSet* normset , const char* range ) const { return m_x.arg().analyticalIntegralWN ( code , normset , range ) ; } -// +// ============================================================================ Int_t Ostap::MoreRooFit::Id::getAnalyticalIntegral ( RooArgSet& allVars , RooArgSet& analVars , const char* range ) const { return m_x.arg().getAnalyticalIntegral ( allVars , analVars , range ) ; } -// +// ============================================================================ Int_t Ostap::MoreRooFit::Id::getAnalyticalIntegralWN ( RooArgSet& allVars , RooArgSet& analVars , @@ -1107,6 +1107,7 @@ Int_t Ostap::MoreRooFit::Id::getAnalyticalIntegralWN { return m_x.arg().getAnalyticalIntegralWN ( allVars , analVars , normset , range ) ; } // ============================================================================ + // ============================================================================ /* constructor from name, title and two pdfs * @param name name @@ -1262,7 +1263,7 @@ Ostap::MoreRooFit::AddDeps::AddDeps ( const AddDeps& right , const char* newname ) : OneVar ( right , newname ) - , m_vlst ( "!vlst" , this , right.m_vlst ) + , m_vlst ( "!vlst" , this , right.m_vlst ) {} // ============================================================================ // destructor @@ -1275,6 +1276,30 @@ Ostap::MoreRooFit::AddDeps* Ostap::MoreRooFit::AddDeps::clone ( const char* newname ) const { return new Ostap::MoreRooFit::AddDeps( *this , newname ) ; } // ============================================================================ +Double_t Ostap::MoreRooFit::AddDeps::analyticalIntegral +( Int_t code , + const char* range ) const +{ return m_x.arg().analyticalIntegral ( code , range ) ; } +// ============================================================================ +Double_t Ostap::MoreRooFit::AddDeps::analyticalIntegralWN +( Int_t code , + const RooArgSet* normset , + const char* range ) const +{ return m_x.arg().analyticalIntegralWN ( code , normset , range ) ; } +// ============================================================================ +Int_t Ostap::MoreRooFit::AddDeps::getAnalyticalIntegral +( RooArgSet& allVars , + RooArgSet& analVars , + const char* range ) const +{ return m_x.arg().getAnalyticalIntegral ( allVars , analVars , range ) ; } +// ============================================================================ +Int_t Ostap::MoreRooFit::AddDeps::getAnalyticalIntegralWN +( RooArgSet& allVars , + RooArgSet& analVars , + const RooArgSet* normset , + const char* range ) const +{ return m_x.arg().getAnalyticalIntegralWN ( allVars , analVars , normset , range ) ; } +// ============================================================================ diff --git a/source/src/MoreVars.cpp b/source/src/MoreVars.cpp index 97226963..450bfed3 100644 --- a/source/src/MoreVars.cpp +++ b/source/src/MoreVars.cpp @@ -44,6 +44,8 @@ namespace const std::string s_v3 = "Ostap::MoreRooFit::Convex" ; const std::string s_v4 = "Ostap::MoreRooFit::ConvexOnly" ; const std::string s_v5 = "Ostap::MoreRooFit::BSpline" ; + const std::string s_v6 = "Ostap::MoreRooFit::Rational" ; + const std::string s_v7 = "Ostap::MoreRooFit::RationalBernstein" ; // =========================================================================== class FakeRecursiveFraction : public RooRecursiveFraction { @@ -273,11 +275,13 @@ namespace // ========================================================================== } // The end of anonymous namespace // ============================================================================ -ClassImp ( Ostap::MoreRooFit::Bernstein ) ; -ClassImp ( Ostap::MoreRooFit::Monotonic ) ; -ClassImp ( Ostap::MoreRooFit::Convex ) ; -ClassImp ( Ostap::MoreRooFit::ConvexOnly ) ; -ClassImp ( Ostap::MoreRooFit::BSpline ) ; +ClassImp ( Ostap::MoreRooFit::Bernstein ) ; +ClassImp ( Ostap::MoreRooFit::Monotonic ) ; +ClassImp ( Ostap::MoreRooFit::Convex ) ; +ClassImp ( Ostap::MoreRooFit::ConvexOnly ) ; +ClassImp ( Ostap::MoreRooFit::BSpline ) ; +ClassImp ( Ostap::MoreRooFit::Rational ) ; +ClassImp ( Ostap::MoreRooFit::RationalBernstein ) ; // ============================================================================ // constructor from the variable, range and list of coefficients // ============================================================================ @@ -357,6 +361,7 @@ Double_t Ostap::MoreRooFit::Bernstein::analyticalIntegral // ============================================================================ + // ============================================================================ // constructor from the variable, range and list of coefficients // ============================================================================ @@ -697,11 +702,194 @@ Double_t Ostap::MoreRooFit::ConvexOnly::analyticalIntegral // return a * ( xmax - xmin ) + b * m_convex.integral ( xmin , xmax ) ; } +// ============================================================================ +// ============================================================================ +// constructor from the variable, range and list of coefficients +// ============================================================================ +Ostap::MoreRooFit::Rational::Rational +( const std::string& name , + const std::string& title , + RooAbsReal& xvar , + const RooArgList& pars , + const unsigned short d , + const double xmin , + const double xmax ) + : RooAbsReal ( name.c_str () , title.c_str () ) + , m_xvar ( "!x" , "Observable" , this , xvar ) + , m_pars ( "!pars" , "Parameters" , this ) + , m_rational ( ::size ( pars ) , d , xmin , xmax ) +{ + // + Ostap::Assert ( 1 <= ::size ( pars ) , s_EMPTYPARS , s_v6 , 510 ) ; + // + ::copy_real ( pars , m_pars , s_INVALIDPAR , s_v6 ) ; + // + Ostap::Assert ( m_rational.npars() == ::size ( m_pars ), s_INVALIDPARS , s_v6 , 512 ) ; + // +} +// ============================================================================= +// copy constructor +// ============================================================================= +Ostap::MoreRooFit::Rational::Rational +( const Ostap::MoreRooFit::Rational& right , + const char* name ) + : RooAbsReal ( right , name ) + , m_xvar ( "!x" , this , right.m_xvar ) + , m_pars ( "!pars" , this , right.m_pars ) + , m_rational ( right.m_rational ) +{} +// ============================================================================= +// default constructor +// ============================================================================ +Ostap::MoreRooFit::Rational::Rational(){} +// ============================================================================= +// destructor +// ============================================================================ +Ostap::MoreRooFit::Rational::~Rational(){} +// ============================================================================ +// clone it! +// ============================================================================ +Ostap::MoreRooFit::Rational* +Ostap::MoreRooFit::Rational::clone ( const char* newname ) const +{ return new Rational( *this , newname ) ; } +// ============================================================================ +void Ostap::MoreRooFit::Rational::setPars() const +{ ::set_pars ( m_pars , m_rational ) ; } +// ============================================================================ +Double_t Ostap::MoreRooFit::Rational::evaluate () const +{ + setPars () ; + const double x = m_xvar ; + return m_rational ( x ) ; +} +// ============================================================================ +Int_t Ostap::MoreRooFit::Rational::getAnalyticalIntegral +( RooArgSet& allVars , + RooArgSet& analVars , + const char* /* rangeName */ ) const +{ + return matchArgs ( allVars , analVars , m_xvar ) ? 1 : 0 ; +} +// =========================================================================== +Double_t Ostap::MoreRooFit::Rational::analyticalIntegral +( Int_t code , + const char* rangeName ) const +{ + setPars() ; + const double xmin = m_xvar.min ( rangeName ) ; + const double xmax = m_xvar.max ( rangeName ) ; + return m_rational.integral ( xmin , xmax ) ; +} +// ============================================================================ + +// ============================================================================ +// constructor from the variable, range and list of coefficients +// ============================================================================ +Ostap::MoreRooFit::RationalBernstein::RationalBernstein +( const std::string& name , + const std::string& title , + RooAbsReal& xvar , + const RooArgList& pars , + const unsigned short p , + const double xmin , + const double xmax ) + : RooAbsReal ( name.c_str () , title.c_str () ) + , m_xvar ( "!x" , "Observable" , this , xvar ) + , m_pars ( "!pars" , "Parameters" , this ) + , m_rational ( ::size ( pars ) , p , xmin , xmax ) +{ + // + Ostap::Assert ( 1 <= ::size ( pars ) , s_EMPTYPARS , s_v7 , 510 ) ; + // + ::copy_real ( pars , m_pars , s_INVALIDPAR , s_v7 ) ; + // + Ostap::Assert ( m_rational.npars() == ::size ( m_pars ), s_INVALIDPARS , s_v7 , 512 ) ; + // +} +// ============================================================================ +// constructor from the variable, range and list of coefficients +// ============================================================================ +Ostap::MoreRooFit::RationalBernstein::RationalBernstein +( const std::string& name , + const std::string& title , + RooAbsReal& xvar , + const RooArgList& p , + const RooArgList& q , + const double xmin , + const double xmax ) + : RooAbsReal ( name.c_str () , title.c_str () ) + , m_xvar ( "!x" , "Observable" , this , xvar ) + , m_pars ( "!pars" , "Parameters" , this ) + , m_rational ( ::size ( p ) , ::size ( q ) , xmin , xmax ) +{ + // + Ostap::Assert ( 1 <= ::size ( p ) + ::size ( q ) , s_EMPTYPARS , s_v7 , 510 ) ; + // + ::copy_real ( p , m_pars , s_INVALIDPAR , s_v7 ) ; + ::copy_real ( q , m_pars , s_INVALIDPAR , s_v7 ) ; + // + Ostap::Assert ( m_rational.npars() == ::size ( m_pars ), s_INVALIDPARS , s_v7 , 512 ) ; + // +} +// ============================================================================= +// copy constructor +// ============================================================================= +Ostap::MoreRooFit::RationalBernstein::RationalBernstein +( const Ostap::MoreRooFit::RationalBernstein& right , + const char* name ) + : RooAbsReal ( right , name ) + , m_xvar ( "!x" , this , right.m_xvar ) + , m_pars ( "!pars" , this , right.m_pars ) + , m_rational ( right.m_rational ) +{} +// ============================================================================= +// default constructor +// ============================================================================ +Ostap::MoreRooFit::RationalBernstein::RationalBernstein(){} +// ============================================================================= +// destructor +// ============================================================================ +Ostap::MoreRooFit::RationalBernstein::~RationalBernstein(){} +// ============================================================================ +// clone it! +// ============================================================================ +Ostap::MoreRooFit::RationalBernstein* +Ostap::MoreRooFit::RationalBernstein::clone ( const char* newname ) const +{ return new RationalBernstein ( *this , newname ) ; } +// ============================================================================ +void Ostap::MoreRooFit::RationalBernstein::setPars() const +{ ::set_pars ( m_pars , m_rational ) ; } +// ============================================================================ +Double_t Ostap::MoreRooFit::RationalBernstein::evaluate () const +{ + setPars () ; + const double x = m_xvar ; + return m_rational ( x ) ; +} +// ============================================================================ +Int_t Ostap::MoreRooFit::RationalBernstein::getAnalyticalIntegral +( RooArgSet& allVars , + RooArgSet& analVars , + const char* /* rangeName */ ) const +{ + return matchArgs ( allVars , analVars , m_xvar ) ? 1 : 0 ; +} +// =========================================================================== +Double_t Ostap::MoreRooFit::RationalBernstein::analyticalIntegral +( Int_t code , + const char* rangeName ) const +{ + setPars() ; + const double xmin = m_xvar.min ( rangeName ) ; + const double xmax = m_xvar.max ( rangeName ) ; + return m_rational.integral ( xmin , xmax ) ; +} +// ============================================================================