diff --git a/roofit/roofitcore/inc/RooFormulaVar.h b/roofit/roofitcore/inc/RooFormulaVar.h index f874b74250fd6..d722153715b6f 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,6 +72,14 @@ class RooFormulaVar : public RooAbsReal { double defaultErrorLevel() const override ; + // Declare this function to be piecewise constant (flat) within the bins of + // the given `binning` of observable `obs`. This lets integration use the fast + // bin integrator instead of the generic numeric integrator. A binning can be + // set for more than one observable. Use a RooUniformBinning to describe many + // uniform bins compactly. + void setBinBoundaries(RooAbsRealLValue &obs, const RooAbsBinning &binning, bool checkFlatness = true); + + 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 ; @@ -94,7 +106,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..39c091829ac26 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,18 @@ class RooGenericPdf : public RooAbsPdf { std::string getUniqueFuncName() const; -protected: + // Declare this pdf to be piecewise constant (flat) within the bins of the + // given `binning` of observable `obs`. This lets the integration use the fast + // bin integrator instead of the generic numeric integrator. A binning can be + // set for more than one observable. Use a RooUniformBinning to describe many + // uniform bins compactly. + void setBinBoundaries(RooAbsRealLValue &obs, const RooAbsBinning &binning, bool checkFlatness = true); + + 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 +94,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..4554d21b20144 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,74 @@ 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::setBinBoundaries() for details. + +void RooFormulaVar::setBinBoundaries(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::setBinBoundaries(" << GetName() << ") the observable " << obs.GetName() + << " is not one of the formula variables of this function, nothing done." << std::endl; + return; + } + + if (checkFlatness) { + // Vary the actual formula variable (the server), which may be a different + // object than `obs` if `obs` is just a same-named stand-in. + auto *serverObs = dynamic_cast(_actualVars.at(idx)); + RooAbsRealLValue &flatObs = serverObs ? *serverObs : obs; + std::span boundaries{binning.array(), static_cast(binning.numBoundaries())}; + if (!RooHelpers::isFunctionFlatInBins(*this, flatObs, boundaries)) { + coutE(InputArguments) << "RooFormulaVar::setBinBoundaries(" << 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()}; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Return true if a binning was set with setBinBoundaries() 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) { + if (_binnings.find(_actualVars.index(o->GetName())) == _binnings.end()) { + return false; + } + } + return true; +} //////////////////////////////////////////////////////////////////////////////// -/// Forward the plot sampling hint from the p.d.f. that defines the observable obs +/// Return the boundaries of the binning set with setBinBoundaries() 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 +315,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..ad1cd39ef1499 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,99 @@ 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::setBinBoundaries(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::setBinBoundaries(" << GetName() << ") the observable " << obs.GetName() + << " is not one of the formula variables of this pdf, nothing done." << std::endl; + return; + } + + if (checkFlatness) { + // Vary the actual formula variable (the server), which may be a different + // object than `obs` if `obs` is just a same-named stand-in. + auto *serverObs = dynamic_cast(_actualVars.at(idx)); + RooAbsRealLValue &flatObs = serverObs ? *serverObs : obs; + std::span boundaries{binning.array(), static_cast(binning.numBoundaries())}; + if (!RooHelpers::isFunctionFlatInBins(*this, flatObs, boundaries)) { + coutE(InputArguments) << "RooGenericPdf::setBinBoundaries(" << 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()}; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Return true if a binning was set with setBinBoundaries() 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) { + if (_binnings.find(_actualVars.index(o->GetName())) == _binnings.end()) { + return false; + } + } + return true; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Return the boundaries of the binning set with setBinBoundaries() 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..e7fcf6df58c18 100644 --- a/roofit/roofitcore/test/testGenericPdf.cxx +++ b/roofit/roofitcore/test/testGenericPdf.cxx @@ -4,11 +4,16 @@ #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 +88,167 @@ TEST(GenericPdf, IdentidyPdfNormalization) EXPECT_NEAR(exp.getVal(), pdf.getVal(), tol); EXPECT_NEAR(exp.getVal(normSet), pdf.getVal(normSet), tol); } + +// Tests for RooGenericPdf::setBinBoundaries(), 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.setBinBoundaries(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); +} + +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.setBinBoundaries(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); +} + +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.setBinBoundaries(x, RooBinning(4, bad)); + EXPECT_FALSE(pdf.isBinnedDistribution(RooArgSet(x))); + + // With the flatness check disabled, the binning is accepted regardless. + pdf.setBinBoundaries(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.setBinBoundaries(other, RooUniformBinning(0.0, 1.0, 5)); + EXPECT_FALSE(pdf.isBinnedDistribution(RooArgSet(other))); +} + +// 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.setBinBoundaries(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.setBinBoundaries(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.setBinBoundaries(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.setBinBoundaries(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); +} diff --git a/roofit/roofitcore/test/testRooFormula.cxx b/roofit/roofitcore/test/testRooFormula.cxx index b250dbc23d953..a57fceb2aabc9 100644 --- a/roofit/roofitcore/test/testRooFormula.cxx +++ b/roofit/roofitcore/test/testRooFormula.cxx @@ -7,11 +7,16 @@ #include #include #include +#include +#include +#include #include #include +#include + /// Since TFormula does very surprising things, /// RooFit needs to do safety checks. /// ``` @@ -109,3 +114,29 @@ TEST(RooFormula, RooConstVarSafeSubstitution) ASSERT_ANY_THROW(RooFormulaVar f1("f1", "x + 0", {x, zero})) << "RooConst variables, if having numeric name, should have name value equal to actual value."; } + +// RooFormulaVar::setBinBoundaries() declares the formula to be piecewise +// constant (flat) within bins of an observable, so that integration uses the +// fast bin integrator instead of the generic numeric integrator. +TEST(RooFormula, SetBinBoundaries) +{ + 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 *formula = "(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"; + RooFormulaVar func("func", "", formula, RooArgList(x, h0, h1, h2, h3, h4)); + + func.setBinBoundaries(x, RooUniformBinning(0.0, 10.0, 5)); + EXPECT_TRUE(func.isBinnedDistribution(RooArgSet(x))); + + std::string integratorName; + double value = 0.0; + { + RooHelpers::HijackMessageStream hijack(RooFit::INFO, RooFit::NumericIntegration); + std::unique_ptr integ{func.createIntegral(x)}; + value = integ->getVal(); + integratorName = hijack.str(); + } + EXPECT_NE(integratorName.find("RooBinIntegrator"), std::string::npos) << integratorName; + EXPECT_DOUBLE_EQ(value, (1.0 + 3.0 + 2.0 + 4.0 + 1.5) * 2.0); // = 23 +}