From 84cccc075e5885a537079ef836bbf6c9cfbfb1b3 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Thu, 25 Jun 2026 19:50:30 +0200 Subject: [PATCH] [RF] Let RooGenericPdf/RooFormulaVar declare a piecewise-flat binning A RooGenericPdf or RooFormulaVar that is constant within bins of an observable (a step function) was always integrated with the generic adaptive numeric integrator, which is needlessly expensive: for a flat distribution the integral is just the sum of each bin's value times its width. Add RooGenericPdf::setBinning() and RooFormulaVar::setBinning(), which take a RooAbsBinning for an observable and declare the function to be flat within those bins. Once set: * isBinnedDistribution() reports the observable as binned, so RooRealIntegral selects the fast RooBinIntegrator; * binBoundaries() and plotSamplingHint() expose the bins so that integration covers exactly the range and plotting draws crisp steps. A binning can be registered for several observables (e.g. for a multi-dimensional flat distribution). By default the function is sampled inside each bin to verify that it really is flat, and the binning is rejected with an error otherwise; pass checkFlatness=false to skip this. The binnings are stored in a std::map keyed by the observable's index in the internal variable list, so they survive renaming of a variable or a server redirection, and a RooUniformBinning keeps the storage compact even for many bins. Lookups resolve the observable by name, so a same-named stand-in (e.g. one read back separately from a file) is accepted too. The shared flatness check and bin-boundary helpers live in RooHelpers. The persistent schema version of both classes is bumped to 2. --- roofit/roofitcore/inc/RooFormulaVar.h | 17 +- roofit/roofitcore/inc/RooGenericPdf.h | 19 +- roofit/roofitcore/inc/RooHelpers.h | 20 ++ roofit/roofitcore/src/RooFormulaVar.cxx | 103 ++++++++- roofit/roofitcore/src/RooGenericPdf.cxx | 119 ++++++++++ roofit/roofitcore/src/RooHelpers.cxx | 58 +++++ roofit/roofitcore/test/testGenericPdf.cxx | 257 ++++++++++++++++++++++ 7 files changed, 584 insertions(+), 9 deletions(-) diff --git a/roofit/roofitcore/inc/RooFormulaVar.h b/roofit/roofitcore/inc/RooFormulaVar.h index f874b74250fd6..da2899f18137d 100644 --- a/roofit/roofitcore/inc/RooFormulaVar.h +++ b/roofit/roofitcore/inc/RooFormulaVar.h @@ -20,12 +20,16 @@ #include "RooArgList.h" #include "RooListProxy.h" #include "RooTrace.h" +#include "RooAbsBinning.h" #include #include +#include +#include class RooArgSet ; class RooFormula ; +class RooAbsRealLValue; class RooFormulaVar : public RooAbsReal { public: @@ -68,8 +72,12 @@ class RooFormulaVar : public RooAbsReal { double defaultErrorLevel() const override ; - std::list* binBoundaries(RooAbsRealLValue& /*obs*/, double /*xlo*/, double /*xhi*/) const override ; - std::list* plotSamplingHint(RooAbsRealLValue& /*obs*/, double /*xlo*/, double /*xhi*/) const override ; + void setBinning(const RooAbsRealLValue &obs, const RooAbsBinning &binning, bool checkFlatness = true); + bool removeBinning(const RooAbsRealLValue &obs); + + bool isBinnedDistribution(const RooArgSet &obs) const override; + std::list* binBoundaries(RooAbsRealLValue& obs, double xlo, double xhi) const override ; + std::list* plotSamplingHint(RooAbsRealLValue& obs, double xlo, double xhi) const override ; // Function evaluation double evaluate() const override ; @@ -94,7 +102,10 @@ class RooFormulaVar : public RooAbsReal { mutable RooArgSet* _nset{nullptr}; ///> _binnings; ///< User-defined binnings, keyed by the observable's index + ///< in _actualVars, for a piecewise-flat distribution + + ClassDefOverride(RooFormulaVar, 2) // Real-valued function of other RooAbsArgs calculated by a TFormula expression }; #endif diff --git a/roofit/roofitcore/inc/RooGenericPdf.h b/roofit/roofitcore/inc/RooGenericPdf.h index e4c779cbd54fc..4270c5c8da5ce 100644 --- a/roofit/roofitcore/inc/RooGenericPdf.h +++ b/roofit/roofitcore/inc/RooGenericPdf.h @@ -18,9 +18,15 @@ #include "RooAbsPdf.h" #include "RooListProxy.h" +#include "RooAbsBinning.h" + +#include +#include +#include class RooArgList ; class RooFormula ; +class RooAbsRealLValue; class RooGenericPdf : public RooAbsPdf { public: @@ -61,8 +67,14 @@ class RooGenericPdf : public RooAbsPdf { std::string getUniqueFuncName() const; -protected: + void setBinning(const RooAbsRealLValue &obs, const RooAbsBinning &binning, bool checkFlatness = true); + bool removeBinning(const RooAbsRealLValue &obs); + + bool isBinnedDistribution(const RooArgSet &obs) const override; + std::list *binBoundaries(RooAbsRealLValue &obs, double xlo, double xhi) const override; + std::list *plotSamplingHint(RooAbsRealLValue &obs, double xlo, double xhi) const override; +protected: RooFormula& formula() const ; // Function evaluation @@ -78,7 +90,10 @@ class RooGenericPdf : public RooAbsPdf { mutable RooFormula * _formula = nullptr; ///> _binnings; ///< User-defined binnings, keyed by the observable's index + ///< in _actualVars, for a piecewise-flat distribution + + ClassDefOverride(RooGenericPdf, 2) // Generic PDF defined by string expression and list of variables }; #endif diff --git a/roofit/roofitcore/inc/RooHelpers.h b/roofit/roofitcore/inc/RooHelpers.h index cb1ebd71c5a11..75362c204beb6 100644 --- a/roofit/roofitcore/inc/RooHelpers.h +++ b/roofit/roofitcore/inc/RooHelpers.h @@ -17,13 +17,17 @@ #include #include +#include + #include +#include #include #include #include class RooAbsPdf; class RooAbsData; +class RooAbsRealLValue; namespace RooHelpers { @@ -91,6 +95,22 @@ void checkRangeOfParameters(const RooAbsReal *callingClass, std::initializer_lis /// set all RooRealVars to constants. return true if at least one changed status bool setAllConstant(const RooAbsCollection &coll, bool constant = true); +/// Check that `function` is constant (flat) inside each bin defined by the +/// sorted `boundaries` when scanning the observable `obs`. Several interior +/// points are sampled per bin and compared to the bin's first sample; if any +/// of them deviates by more than `relTol` (relative to the value scale), the +/// function is not flat and false is returned. The value of `obs` is restored +/// on return. +bool isFunctionFlatInBins(const RooAbsReal &function, RooAbsRealLValue &obs, std::span boundaries, + double relTol = 1e-9); + +/// Return a newly allocated list with the subset of `boundaries` that lies +/// strictly inside [`xlo`, `xhi`], with `xlo` and `xhi` added as the first and +/// last entries. This is the form expected by RooFit's binBoundaries() +/// interface, so the bin integrator covers exactly the integration range. +/// The caller takes ownership of the returned list. +std::list *binBoundariesInRange(std::span boundaries, double xlo, double xhi); + } // namespace RooHelpers #endif diff --git a/roofit/roofitcore/src/RooFormulaVar.cxx b/roofit/roofitcore/src/RooFormulaVar.cxx index c6fd11cbb7b4b..b1e19db40f145 100644 --- a/roofit/roofitcore/src/RooFormulaVar.cxx +++ b/roofit/roofitcore/src/RooFormulaVar.cxx @@ -52,6 +52,10 @@ #include "RooMsgService.h" #include "RooTrace.h" #include "RooFormula.h" +#include "RooAbsRealLValue.h" +#include "RooAbsBinning.h" +#include "RooCurve.h" +#include "RooHelpers.h" #ifdef ROOFIT_LEGACY_EVAL_BACKEND #include "RooNLLVar.h" @@ -123,6 +127,9 @@ RooFormulaVar::RooFormulaVar(const RooFormulaVar& other, const char* name) : _actualVars("actualVars",this,other._actualVars), _formExpr(other._formExpr) { + for (auto const &item : other._binnings) { + _binnings[item.first] = std::unique_ptr{item.second->clone()}; + } if (other._formula && other._formula->ok()) { _formula = new RooFormula(*other._formula); _formExpr = _formula->formulaString().c_str(); @@ -228,13 +235,94 @@ void RooFormulaVar::writeToStream(ostream& os, bool compact) const } } +//////////////////////////////////////////////////////////////////////////////// +/// Declare that this function is piecewise constant (flat) within the bins of +/// the given `binning` of the observable `obs`, which must be one of the formula +/// variables. The method can be called several times to set a binning for more +/// than one observable. See RooGenericPdf::setBinning() for details. +void RooFormulaVar::setBinning(const RooAbsRealLValue &obs, const RooAbsBinning &binning, bool checkFlatness) +{ + // Match the observable to a formula variable by name, so that a same-named + // stand-in for the actual server is accepted too. + const int idx = _actualVars.index(obs.GetName()); + if (idx < 0) { + coutE(InputArguments) << "RooFormulaVar::setBinning(" << GetName() << ") the observable " << obs.GetName() + << " is not one of the formula variables of this function, nothing done." << std::endl; + return; + } + + if (checkFlatness) { + // Sample the function by varying the actual formula variable (the server), + // which may be a different object than `obs` if `obs` is just a same-named + // stand-in: the function's value depends on the server, not on `obs`. + if (auto *serverObs = dynamic_cast(_actualVars.at(idx))) { + std::span boundaries{binning.array(), static_cast(binning.numBoundaries())}; + if (!RooHelpers::isFunctionFlatInBins(*this, *serverObs, boundaries)) { + coutE(InputArguments) << "RooFormulaVar::setBinning(" << GetName() << ") the expression \"" << _formExpr + << "\" is not flat within the given bins of " << obs.GetName() + << ". The binning is not set. Pass checkFlatness=false to override this check." + << std::endl; + return; + } + } + } + + // Key the binning by the observable's index in _actualVars (not its name), so + // that it survives a renaming of the variable or a server redirection. + _binnings[idx] = std::unique_ptr{binning.clone()}; +} //////////////////////////////////////////////////////////////////////////////// -/// Forward the plot sampling hint from the p.d.f. that defines the observable obs +/// Remove a binning previously declared with setBinning() for observable `obs`, +/// reverting to the generic numeric integrator for it. Returns true if a binning +/// was removed, false if none was set for `obs`. + +bool RooFormulaVar::removeBinning(const RooAbsRealLValue &obs) +{ + return _binnings.erase(_actualVars.index(obs.GetName())) > 0; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Return true if a binning was set with setBinning() for every +/// observable in the integration set `obs`. + +bool RooFormulaVar::isBinnedDistribution(const RooArgSet &obs) const +{ + if (obs.empty() || _binnings.empty()) { + return false; + } + for (RooAbsArg *o : obs) { + const int idx = _actualVars.index(o->GetName()); + // Observables that are not formula variables of this function are ones we + // do not depend on: the function is constant (hence trivially binned) in + // them, so they must be ignored here. This matches the convention that + // composite functions like RooProduct rely on, where each component's + // isBinnedDistribution() is queried with the full observable set. + if (idx < 0) { + continue; + } + if (_binnings.find(idx) == _binnings.end()) { + return false; + } + } + return true; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Return the boundaries of the binning set with setBinning() that fall +/// within [xlo, xhi]. If no binning was set for this observable, forward the bin +/// boundaries from the server that defines the observable obs. std::list* RooFormulaVar::binBoundaries(RooAbsRealLValue& obs, double xlo, double xhi) const { + auto found = _binnings.find(_actualVars.index(obs.GetName())); + if (found != _binnings.end()) { + const RooAbsBinning &binning = *found->second; + return RooHelpers::binBoundariesInRange({binning.array(), static_cast(binning.numBoundaries())}, xlo, + xhi); + } + for (const auto par : _actualVars) { auto func = static_cast(par); list* binb = nullptr; @@ -247,13 +335,20 @@ std::list* RooFormulaVar::binBoundaries(RooAbsRealLValue& obs, double xl return nullptr; } - - //////////////////////////////////////////////////////////////////////////////// -/// Forward the plot sampling hint from the p.d.f. that defines the observable obs +/// Return sampling hints that draw the piecewise-flat shape exactly if a binning +/// was set for this observable. Otherwise, forward the plot sampling hint from +/// the server that defines the observable obs. std::list* RooFormulaVar::plotSamplingHint(RooAbsRealLValue& obs, double xlo, double xhi) const { + auto found = _binnings.find(_actualVars.index(obs.GetName())); + if (found != _binnings.end()) { + const RooAbsBinning &binning = *found->second; + return RooCurve::plotSamplingHintForBinBoundaries( + {binning.array(), static_cast(binning.numBoundaries())}, xlo, xhi); + } + for (const auto par : _actualVars) { auto func = dynamic_cast(par); list* hint = nullptr; diff --git a/roofit/roofitcore/src/RooGenericPdf.cxx b/roofit/roofitcore/src/RooGenericPdf.cxx index 3a87554d35ee7..874054865dfab 100644 --- a/roofit/roofitcore/src/RooGenericPdf.cxx +++ b/roofit/roofitcore/src/RooGenericPdf.cxx @@ -48,6 +48,10 @@ the names of the arguments are not hard coded. #include "RooMsgService.h" #include "RooArgList.h" #include "RooFormula.h" +#include "RooAbsRealLValue.h" +#include "RooAbsBinning.h" +#include "RooCurve.h" +#include "RooHelpers.h" using std::istream, std::ostream, std::endl; @@ -107,6 +111,9 @@ RooGenericPdf::RooGenericPdf(const RooGenericPdf& other, const char* name) : _actualVars("actualVars",this,other._actualVars), _formExpr(other._formExpr) { + for (auto const &item : other._binnings) { + _binnings[item.first] = std::unique_ptr{item.second->clone()}; + } formula(); } @@ -122,7 +129,119 @@ RooFormula& RooGenericPdf::formula() const return *_formula ; } +//////////////////////////////////////////////////////////////////////////////// +/// Declare that this pdf is piecewise constant (flat) within the bins of the +/// given `binning` of the observable `obs`, which must be one of the formula +/// variables. The method can be called several times to set a binning for more +/// than one observable. Use a RooUniformBinning to describe many uniform bins +/// compactly. +/// +/// Once set, integrals over `obs` use the fast bin integrator (which sums the +/// central value of each bin times the bin width) instead of the generic +/// numeric integrator, and plotting samples the step shape exactly. +/// +/// If `checkFlatness` is true (the default), the function is sampled at several +/// points inside each bin to verify that it is indeed flat; if it is not, an +/// error is issued and the binning is not stored. + +void RooGenericPdf::setBinning(const RooAbsRealLValue &obs, const RooAbsBinning &binning, bool checkFlatness) +{ + // Match the observable to a formula variable by name, so that a same-named + // stand-in for the actual server is accepted too. + const int idx = _actualVars.index(obs.GetName()); + if (idx < 0) { + coutE(InputArguments) << "RooGenericPdf::setBinning(" << GetName() << ") the observable " << obs.GetName() + << " is not one of the formula variables of this pdf, nothing done." << std::endl; + return; + } + + if (checkFlatness) { + // Sample the function by varying the actual formula variable (the server), + // which may be a different object than `obs` if `obs` is just a same-named + // stand-in: the function's value depends on the server, not on `obs`. + if (auto *serverObs = dynamic_cast(_actualVars.at(idx))) { + std::span boundaries{binning.array(), static_cast(binning.numBoundaries())}; + if (!RooHelpers::isFunctionFlatInBins(*this, *serverObs, boundaries)) { + coutE(InputArguments) << "RooGenericPdf::setBinning(" << GetName() << ") the expression \"" << _formExpr + << "\" is not flat within the given bins of " << obs.GetName() + << ". The binning is not set. Pass checkFlatness=false to override this check." + << std::endl; + return; + } + } + } + + // Key the binning by the observable's index in _actualVars (not its name), so + // that it survives a renaming of the variable or a server redirection. + _binnings[idx] = std::unique_ptr{binning.clone()}; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Remove a binning previously declared with setBinning() for observable `obs`, +/// reverting to the generic numeric integrator for it. Returns true if a binning +/// was removed, false if none was set for `obs`. + +bool RooGenericPdf::removeBinning(const RooAbsRealLValue &obs) +{ + return _binnings.erase(_actualVars.index(obs.GetName())) > 0; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Return true if a binning was set with setBinning() for every +/// observable in the integration set `obs`. + +bool RooGenericPdf::isBinnedDistribution(const RooArgSet &obs) const +{ + if (obs.empty() || _binnings.empty()) { + return false; + } + for (RooAbsArg *o : obs) { + const int idx = _actualVars.index(o->GetName()); + // Observables that are not formula variables of this pdf are ones we do + // not depend on: the function is constant (hence trivially binned) in + // them, so they must be ignored here. This matches the convention that + // composite functions like RooProduct rely on, where each component's + // isBinnedDistribution() is queried with the full observable set. + if (idx < 0) { + continue; + } + if (_binnings.find(idx) == _binnings.end()) { + return false; + } + } + return true; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Return the boundaries of the binning set with setBinning() that fall +/// within [xlo, xhi], or a null pointer if no binning was set for this observable. +std::list *RooGenericPdf::binBoundaries(RooAbsRealLValue &obs, double xlo, double xhi) const +{ + auto found = _binnings.find(_actualVars.index(obs.GetName())); + if (found == _binnings.end()) { + return nullptr; + } + const RooAbsBinning &binning = *found->second; + return RooHelpers::binBoundariesInRange({binning.array(), static_cast(binning.numBoundaries())}, xlo, + xhi); +} + +//////////////////////////////////////////////////////////////////////////////// +/// Return sampling hints that draw the piecewise-flat shape exactly (a pair of +/// points just left and right of every bin boundary), or a null pointer if no +/// binning was set for this observable. + +std::list *RooGenericPdf::plotSamplingHint(RooAbsRealLValue &obs, double xlo, double xhi) const +{ + auto found = _binnings.find(_actualVars.index(obs.GetName())); + if (found == _binnings.end()) { + return nullptr; + } + const RooAbsBinning &binning = *found->second; + return RooCurve::plotSamplingHintForBinBoundaries( + {binning.array(), static_cast(binning.numBoundaries())}, xlo, xhi); +} //////////////////////////////////////////////////////////////////////////////// /// Calculate current value of this object diff --git a/roofit/roofitcore/src/RooHelpers.cxx b/roofit/roofitcore/src/RooHelpers.cxx index 322a451edc137..766170309ca4b 100644 --- a/roofit/roofitcore/src/RooHelpers.cxx +++ b/roofit/roofitcore/src/RooHelpers.cxx @@ -29,6 +29,8 @@ #include #include +#include +#include #include namespace RooHelpers { @@ -168,4 +170,60 @@ bool setAllConstant(const RooAbsCollection &coll, bool constant) return changed; } +bool isFunctionFlatInBins(const RooAbsReal &function, RooAbsRealLValue &obs, std::span boundaries, + double relTol) +{ + // Fractions of the bin width at which the function is sampled. They are kept + // strictly inside the bin (away from the boundaries) so that the evaluation + // is not affected by which side of a boundary a step function jumps. + const double fractions[] = {0.04, 0.27, 0.5, 0.73, 0.96}; + + const double savedVal = obs.getVal(); + + bool isFlat = true; + for (std::size_t i = 0; i + 1 < boundaries.size() && isFlat; ++i) { + const double lo = boundaries[i]; + const double hi = boundaries[i + 1]; + double reference = 0.0; + bool first = true; + for (double frac : fractions) { + obs.setVal(lo + frac * (hi - lo)); + const double val = function.getVal(); + if (first) { + reference = val; + first = false; + continue; + } + const double scale = std::max(std::abs(reference), 1e-12); + if (std::abs(val - reference) > relTol * scale) { + isFlat = false; + break; + } + } + } + + obs.setVal(savedVal); + return isFlat; +} + +std::list *binBoundariesInRange(std::span boundaries, double xlo, double xhi) +{ + auto out = new std::list; + + // Small tolerance so that boundaries numerically coinciding with the range + // limits are not duplicated by the explicit xlo/xhi endpoints below. + const double delta = (xhi - xlo) * 1e-8; + + for (double boundary : boundaries) { + if (boundary > xlo + delta && boundary < xhi - delta) { + out->push_back(boundary); + } + } + + out->push_front(xlo); + out->push_back(xhi); + + return out; +} + } // namespace RooHelpers diff --git a/roofit/roofitcore/test/testGenericPdf.cxx b/roofit/roofitcore/test/testGenericPdf.cxx index 6b5008577dd19..1a3ce3275932e 100644 --- a/roofit/roofitcore/test/testGenericPdf.cxx +++ b/roofit/roofitcore/test/testGenericPdf.cxx @@ -4,11 +4,17 @@ #include #include +#include +#include +#include +#include #include #include #include +#include + #define MAKE_JOHNSON_AND_VARS RooRealVar mass("mass", "mass", 0., -200., 200.);\ RooRealVar mu("mu", "Location parameter of normal distribution", 100., -200., 200.);\ RooRealVar sigma("sigma", "Two sigma of normal distribution", 2., 0., 100.);\ @@ -83,3 +89,254 @@ TEST(GenericPdf, IdentidyPdfNormalization) EXPECT_NEAR(exp.getVal(), pdf.getVal(), tol); EXPECT_NEAR(exp.getVal(normSet), pdf.getVal(normSet), tol); } + +// Tests for RooGenericPdf::setBinning(), which declares the formula to be +// piecewise constant (flat) within bins of an observable so that integration +// can use the fast bin integrator instead of the generic numeric integrator. + +// A pdf that is flat within five uniform bins on [0, 10], with one height per bin. +#define MAKE_PIECEWISE_FLAT_VARS \ + RooRealVar x("x", "x", 0., 10.); \ + RooRealVar h0("h0", "", 1.0), h1("h1", "", 3.0), h2("h2", "", 2.0), h3("h3", "", 4.0), h4("h4", "", 1.5); + +const char *piecewiseFlatFormula = "(floor(x/2)==0)*h0+(floor(x/2)==1)*h1+(floor(x/2)==2)*h2" + "+(floor(x/2)==3)*h3+(floor(x/2)==4)*h4"; + +// Exact integral of the piecewise-flat function: sum(height) * binWidth. +constexpr double piecewiseFlatIntegral = (1.0 + 3.0 + 2.0 + 4.0 + 1.5) * 2.0; // = 23 + +namespace { + +// Integrate `pdf` over `obs`, returning the integral value and reporting the +// name of the numeric integrator that was used via `integratorName`. +double integrate(RooAbsReal &pdf, const RooArgSet &obs, std::string &integratorName) +{ + RooHelpers::HijackMessageStream hijack(RooFit::INFO, RooFit::NumericIntegration); + std::unique_ptr integ{pdf.createIntegral(obs)}; + const double value = integ->getVal(); + integratorName = hijack.str(); + return value; +} + +} // namespace + +TEST(GenericPdf, BinnedIntegrationExplicitBoundaries) +{ + MAKE_PIECEWISE_FLAT_VARS + RooGenericPdf pdf("pdf", "", piecewiseFlatFormula, RooArgList(x, h0, h1, h2, h3, h4)); + + // Without a binning, the generic numeric integrator is used. + { + std::string integratorName; + integrate(pdf, x, integratorName); + EXPECT_EQ(integratorName.find("RooBinIntegrator"), std::string::npos) << integratorName; + } + + const double boundaries[] = {0, 2, 4, 6, 8, 10}; + pdf.setBinning(x, RooBinning(5, boundaries)); + + EXPECT_TRUE(pdf.isBinnedDistribution(RooArgSet(x))); + + std::string integratorName; + const double value = integrate(pdf, x, integratorName); + EXPECT_NE(integratorName.find("RooBinIntegrator"), std::string::npos) << integratorName; + EXPECT_DOUBLE_EQ(value, piecewiseFlatIntegral); +} + +// removeBinning() reverses setBinning(), reverting to the generic integrator. +TEST(GenericPdf, BinnedRemoveBinning) +{ + MAKE_PIECEWISE_FLAT_VARS + RooGenericPdf pdf("pdf", "", piecewiseFlatFormula, RooArgList(x, h0, h1, h2, h3, h4)); + + pdf.setBinning(x, RooUniformBinning(0.0, 10.0, 5)); + ASSERT_TRUE(pdf.isBinnedDistribution(RooArgSet(x))); + + // Removing the binning returns true and reverts to the generic integrator. + EXPECT_TRUE(pdf.removeBinning(x)); + EXPECT_FALSE(pdf.isBinnedDistribution(RooArgSet(x))); + + std::string integratorName; + const double value = integrate(pdf, x, integratorName); + EXPECT_EQ(integratorName.find("RooBinIntegrator"), std::string::npos) << integratorName; + // Same integral, but now via the (approximate) generic numeric integrator. + EXPECT_NEAR(value, piecewiseFlatIntegral, 1e-2); + + // Removing again (or for an unset/unknown observable) returns false. + EXPECT_FALSE(pdf.removeBinning(x)); + RooRealVar other("other", "", 0, 1); + EXPECT_FALSE(pdf.removeBinning(other)); +} + +TEST(GenericPdf, BinnedIntegrationUniformBinning) +{ + MAKE_PIECEWISE_FLAT_VARS + RooGenericPdf pdf("pdf", "", piecewiseFlatFormula, RooArgList(x, h0, h1, h2, h3, h4)); + + // A RooUniformBinning describes the same bins compactly (just min, max, n). + pdf.setBinning(x, RooUniformBinning(0.0, 10.0, 5)); + + EXPECT_TRUE(pdf.isBinnedDistribution(RooArgSet(x))); + + std::string integratorName; + const double value = integrate(pdf, x, integratorName); + EXPECT_NE(integratorName.find("RooBinIntegrator"), std::string::npos) << integratorName; + EXPECT_DOUBLE_EQ(value, piecewiseFlatIntegral); +} + +// The existing tests above all use uniform bins. This one uses bins of +// unequal width to check that a non-uniform binning is integrated correctly: +// the integral must weight each bin's height by its own width. +TEST(GenericPdf, BinnedIntegrationNonUniformBinning) +{ + RooRealVar x("x", "x", 0., 10.); + RooRealVar b0("b0", "", 1.0), b1("b1", "", 3.0), b2("b2", "", 2.0), b3("b3", "", 4.0); + + // A step function that is flat within the non-uniform bins + // [0,1), [1,3), [3,6), [6,10] (widths 1, 2, 3, 4). + const char *formula = "(x<1)*b0 + (x>=1)*(x<3)*b1 + (x>=3)*(x<6)*b2 + (x>=6)*b3"; + RooGenericPdf pdf("pdf", "", formula, RooArgList(x, b0, b1, b2, b3)); + + const double boundaries[] = {0, 1, 3, 6, 10}; + pdf.setBinning(x, RooBinning(4, boundaries)); + + EXPECT_TRUE(pdf.isBinnedDistribution(RooArgSet(x))); + + std::string integratorName; + const double value = integrate(pdf, x, integratorName); + EXPECT_NE(integratorName.find("RooBinIntegrator"), std::string::npos) << integratorName; + // sum(height * width) = 1*1 + 3*2 + 2*3 + 4*4 = 29. + EXPECT_DOUBLE_EQ(value, 1.0 * 1.0 + 3.0 * 2.0 + 2.0 * 3.0 + 4.0 * 4.0); +} + +TEST(GenericPdf, BinnedFlatnessCheck) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::FATAL); // silence the expected error + MAKE_PIECEWISE_FLAT_VARS + RooGenericPdf pdf("pdf", "", piecewiseFlatFormula, RooArgList(x, h0, h1, h2, h3, h4)); + + // Bins that straddle a jump ([1,3], [3,5], ... cross the steps at 2, 4, ...) + // are not flat: the binning must be rejected. + const double bad[] = {1, 3, 5, 7, 9}; + pdf.setBinning(x, RooBinning(4, bad)); + EXPECT_FALSE(pdf.isBinnedDistribution(RooArgSet(x))); + + // With the flatness check disabled, the binning is accepted regardless. + pdf.setBinning(x, RooBinning(4, bad), /*checkFlatness=*/false); + EXPECT_TRUE(pdf.isBinnedDistribution(RooArgSet(x))); +} + +TEST(GenericPdf, BinnedNonFormulaVariable) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::FATAL); // silence the expected error + MAKE_PIECEWISE_FLAT_VARS + RooGenericPdf pdf("pdf", "", piecewiseFlatFormula, RooArgList(x, h0, h1, h2, h3, h4)); + + // An observable that is not one of the formula variables must be rejected. + RooRealVar other("other", "", 0, 1); + pdf.setBinning(other, RooUniformBinning(0.0, 1.0, 5)); + EXPECT_FALSE(pdf.isBinnedDistribution(RooArgSet(other))); +} + +// isBinnedDistribution() must ignore observables that the pdf does not depend +// on: the function is constant, hence trivially binned, in them. Otherwise an +// unrelated observable in the query set would spuriously report the pdf as +// unbinned. +TEST(GenericPdf, BinnedIgnoresUnrelatedObservable) +{ + MAKE_PIECEWISE_FLAT_VARS + RooGenericPdf pdf("pdf", "", piecewiseFlatFormula, RooArgList(x, h0, h1, h2, h3, h4)); + pdf.setBinning(x, RooUniformBinning(0.0, 10.0, 5)); + + RooRealVar y("y", "y", 0., 4.); // not a formula variable of the pdf + EXPECT_TRUE(pdf.isBinnedDistribution(RooArgSet(x, y))); +} + +// A RooGenericPdf that is flat-binned in its own observable must work as a +// binned component of a RooProduct. RooProduct queries each component with the +// full observable set, so each component must ignore the observable it does not +// depend on (see BinnedIgnoresUnrelatedObservable) for the product to report as +// binned. +TEST(GenericPdf, BinnedProductComposition) +{ + RooRealVar x("x", "x", 0., 10.); + RooRealVar y("y", "y", 0., 4.); + RooRealVar a("a", "", 1.0), b("b", "", 2.0); + + RooGenericPdf px("px", "", "a*(x<100)", RooArgList(x, a)); + RooGenericPdf py("py", "", "b*(y<100)", RooArgList(y, b)); + px.setBinning(x, RooUniformBinning(0.0, 10.0, 5)); + py.setBinning(y, RooUniformBinning(0.0, 4.0, 2)); + + EXPECT_TRUE(px.isBinnedDistribution(RooArgSet(x, y))); + EXPECT_TRUE(py.isBinnedDistribution(RooArgSet(x, y))); + + RooProduct prod("prod", "", RooArgList(px, py)); + EXPECT_TRUE(prod.isBinnedDistribution(RooArgSet(x, y))); +} + +// The lookup resolves the observable by name, so a different object with the +// same name (e.g. one read back separately from a file) resolves to the same +// binning - it need not be the pdf's own internal server object. +TEST(GenericPdf, BinnedLookupBySameNamedVariable) +{ + MAKE_PIECEWISE_FLAT_VARS + RooGenericPdf pdf("pdf", "", piecewiseFlatFormula, RooArgList(x, h0, h1, h2, h3, h4)); + pdf.setBinning(x, RooUniformBinning(0.0, 10.0, 5)); + + RooRealVar xStandIn("x", "x", 0., 10.); + EXPECT_NE(&xStandIn, pdf.getParameter("x")); + + EXPECT_TRUE(pdf.isBinnedDistribution(RooArgSet(xStandIn))); + + std::unique_ptr> boundaries{pdf.binBoundaries(xStandIn, 0.0, 10.0)}; + ASSERT_NE(boundaries, nullptr); + EXPECT_EQ(boundaries->size(), 6u); // 5 bins -> 6 boundaries +} + +TEST(GenericPdf, BinnedMultipleObservables) +{ + MAKE_PIECEWISE_FLAT_VARS + // Second observable, flat in two bins on [0, 4]. + RooRealVar y("y", "y", 0., 4.); + RooRealVar g0("g0", "", 2.0), g1("g1", "", 5.0); + + TString formula = TString::Format("(%s) * ((floor(y/2)==0)*g0 + (floor(y/2)==1)*g1)", piecewiseFlatFormula); + RooGenericPdf pdf("pdf", "", formula, RooArgList(x, y, h0, h1, h2, h3, h4, g0, g1)); + + pdf.setBinning(x, RooUniformBinning(0.0, 10.0, 5)); + + // Only x is binned so far: the joint {x,y} distribution is not fully binned. + EXPECT_FALSE(pdf.isBinnedDistribution(RooArgSet(x, y))); + EXPECT_TRUE(pdf.isBinnedDistribution(RooArgSet(x))); + + pdf.setBinning(y, RooUniformBinning(0.0, 4.0, 2)); + EXPECT_TRUE(pdf.isBinnedDistribution(RooArgSet(x, y))); + + std::string integratorName; + const double value = integrate(pdf, RooArgSet(x, y), integratorName); + EXPECT_NE(integratorName.find("RooBinIntegrator"), std::string::npos) << integratorName; + // int_x g(x) dx = 23, int_y g(y) dy = (2+5)*2 = 14, product = 322. + EXPECT_DOUBLE_EQ(value, piecewiseFlatIntegral * ((2.0 + 5.0) * 2.0)); +} + +// The binning is keyed by the observable's index in the internal variable list, +// not by its name, so renaming the observable after registering the binning +// must not lose it. +TEST(GenericPdf, BinnedBoundariesSurviveRename) +{ + MAKE_PIECEWISE_FLAT_VARS + RooGenericPdf pdf("pdf", "", piecewiseFlatFormula, RooArgList(x, h0, h1, h2, h3, h4)); + pdf.setBinning(x, RooUniformBinning(0.0, 10.0, 5)); + ASSERT_TRUE(pdf.isBinnedDistribution(RooArgSet(x))); + + // Rename the observable in place. + x.SetName("xRenamed"); + + EXPECT_TRUE(pdf.isBinnedDistribution(RooArgSet(x))); + + std::string integratorName; + const double value = integrate(pdf, RooArgSet(x), integratorName); + EXPECT_NE(integratorName.find("RooBinIntegrator"), std::string::npos) << integratorName; + EXPECT_DOUBLE_EQ(value, piecewiseFlatIntegral); +}