From 3b38fed2214eeca1a7028f0c9fd1def629345281 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Wed, 20 Nov 2024 18:22:53 +0100 Subject: [PATCH 1/2] [RF] Don't randomize constant parameters in testRooFuncWrapper --- roofit/roofitcore/test/testRooFuncWrapper.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/roofit/roofitcore/test/testRooFuncWrapper.cxx b/roofit/roofitcore/test/testRooFuncWrapper.cxx index 6a372b7daae8f..1146975e3aa5a 100644 --- a/roofit/roofitcore/test/testRooFuncWrapper.cxx +++ b/roofit/roofitcore/test/testRooFuncWrapper.cxx @@ -76,7 +76,7 @@ void randomizeParameters(const RooArgSet ¶meters) double mul = unif(re); auto par = dynamic_cast(param); - if (!par) + if (!par || par->isConstant()) continue; double val = par->getVal(); val = val + mul * (mul > 0 ? (par->getMax() - val) : (val - par->getMin())); @@ -241,6 +241,7 @@ TEST_P(FactoryTest, NLLFit) // We don't use the RooFit::Evaluator for the nominal likelihood. Like this, // we make sure to validate also the NLL values of the generated code. static_cast(*nllFunc).disableEvaluator(); + static_cast(*nllFunc).writeDebugMacro(_params._name); double tol = _params._fitResultTolerance; From d5e95ac1c1c67973157db5704b757a903d878ac1 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Wed, 20 Nov 2024 18:26:38 +0100 Subject: [PATCH 2/2] [RF] Implement channel masking for simultaneous likelihoods Upstreaming a feature from CMS combine. Draft for now. --- roofit/codegen/inc/RooFit/CodegenImpl.h | 2 + roofit/codegen/src/CodegenImpl.cxx | 32 +++++++--- roofit/roofitcore/inc/LinkDef.h | 1 + .../inc/RooFit/Detail/RooNLLVarNew.h | 27 +++++++++ roofit/roofitcore/src/FitHelpers.cxx | 3 +- roofit/roofitcore/src/RooNLLVarNew.cxx | 59 ++++++++++++++++++- 6 files changed, 111 insertions(+), 13 deletions(-) diff --git a/roofit/codegen/inc/RooFit/CodegenImpl.h b/roofit/codegen/inc/RooFit/CodegenImpl.h index a0d22110a13e1..45bdca31e23e3 100644 --- a/roofit/codegen/inc/RooFit/CodegenImpl.h +++ b/roofit/codegen/inc/RooFit/CodegenImpl.h @@ -67,6 +67,7 @@ namespace Detail { class RooFixedProdPdf; class RooNLLVarNew; class RooNormalizedPdf; +class RooSimNLL; } // namespace Detail namespace Experimental { @@ -76,6 +77,7 @@ class CodegenContext; void codegenImpl(RooFit::Detail::RooFixedProdPdf &arg, CodegenContext &ctx); void codegenImpl(RooFit::Detail::RooNLLVarNew &arg, CodegenContext &ctx); void codegenImpl(RooFit::Detail::RooNormalizedPdf &arg, CodegenContext &ctx); +void codegenImpl(RooFit::Detail::RooSimNLL &arg, CodegenContext &ctx); void codegenImpl(ParamHistFunc &arg, CodegenContext &ctx); void codegenImpl(PiecewiseInterpolation &arg, CodegenContext &ctx); void codegenImpl(RooAbsArg &arg, CodegenContext &ctx); diff --git a/roofit/codegen/src/CodegenImpl.cxx b/roofit/codegen/src/CodegenImpl.cxx index a4ef24b67f499..d761a20ae9776 100644 --- a/roofit/codegen/src/CodegenImpl.cxx +++ b/roofit/codegen/src/CodegenImpl.cxx @@ -133,6 +133,28 @@ void codegenImpl(RooFit::Detail::RooFixedProdPdf &arg, CodegenContext &ctx) } } +void codegenImpl(RooFit::Detail::RooSimNLL &arg, CodegenContext &ctx) +{ + if (arg.terms().empty()) { + ctx.addResult(&arg, "0.0"); + } + + std::string resName = RooFit::Detail::makeValidVarName(arg.GetName()) + "Result"; + ctx.addResult(&arg, resName); + ctx.addToGlobalScope("double " + resName + " = 0.0;\n"); + + std::stringstream ss; + + std::size_t i = 0; + for (auto *component : static_range_cast(arg.terms())) { + + // TODO: support channel masking here + ss << resName << " += " << ctx.buildFunction(*component, ctx.outputSizes()) << "(params, obs, xlArr);\n"; + ++i; + } + ctx.addToGlobalScope(ss.str()); +} + void codegenImpl(ParamHistFunc &arg, CodegenContext &ctx) { std::string const &idx = arg.dataHist().calculateTreeIndexForCodeSquash(&arg, ctx, arg.dataVars(), true); @@ -251,15 +273,7 @@ void codegenImpl(RooAddition &arg, CodegenContext &ctx) std::size_t i = 0; for (auto *component : static_range_cast(arg.list())) { - - if (!dynamic_cast(component) || arg.list().size() == 1) { - result += ctx.getResult(*component); - ++i; - if (i < arg.list().size()) - result += '+'; - continue; - } - result += ctx.buildFunction(*component, ctx.outputSizes()) + "(params, obs, xlArr)"; + result += ctx.getResult(*component); ++i; if (i < arg.list().size()) result += '+'; diff --git a/roofit/roofitcore/inc/LinkDef.h b/roofit/roofitcore/inc/LinkDef.h index 83fa0c52e8ac5..9f257896e1ad7 100644 --- a/roofit/roofitcore/inc/LinkDef.h +++ b/roofit/roofitcore/inc/LinkDef.h @@ -337,5 +337,6 @@ #pragma link C++ class RooBinWidthFunction+; #pragma link C++ class RooFit::Detail::RooNLLVarNew+; #pragma link C++ class RooFit::Detail::RooNormalizedPdf+ ; +#pragma link C++ class RooFit::Detail::RooSimNLL+; #endif diff --git a/roofit/roofitcore/inc/RooFit/Detail/RooNLLVarNew.h b/roofit/roofitcore/inc/RooFit/Detail/RooNLLVarNew.h index d7ef8e8fda613..678a00c22f610 100644 --- a/roofit/roofitcore/inc/RooFit/Detail/RooNLLVarNew.h +++ b/roofit/roofitcore/inc/RooFit/Detail/RooNLLVarNew.h @@ -19,10 +19,13 @@ #include #include #include +#include #include #include +class RooAbsCategory; + namespace RooFit { namespace Detail { @@ -87,6 +90,30 @@ class RooNLLVarNew : public RooAbsReal { ClassDefOverride(RooFit::Detail::RooNLLVarNew, 0); }; +class RooSimNLL : public RooAbsReal { +public: + RooSimNLL(const char *name, const char *title, const RooArgSet &terms, RooAbsCategoryLValue const &indexCat, + bool channelMasking); + + RooSimNLL(const RooSimNLL &other, const char *name = nullptr); + TObject *clone(const char *newname) const override { return new RooSimNLL(*this, newname); } + + double defaultErrorLevel() const override; + + const RooArgSet &terms() const { return _set; } + const RooArgSet &masks() const { return _mask; } + + void doEval(RooFit::EvalContext &) const override; + +protected: + double evaluate() const override; + + RooSetProxy _set; ///< set of terms to be summed + RooSetProxy _mask; + + ClassDefOverride(RooFit::Detail::RooSimNLL, 0) // Sum of RooNLLVarNew instances +}; + } // namespace Detail } // namespace RooFit diff --git a/roofit/roofitcore/src/FitHelpers.cxx b/roofit/roofitcore/src/FitHelpers.cxx index 8e5b45ac0c86f..c67d157df852c 100644 --- a/roofit/roofitcore/src/FitHelpers.cxx +++ b/roofit/roofitcore/src/FitHelpers.cxx @@ -51,6 +51,7 @@ #endif using RooFit::Detail::RooNLLVarNew; +using RooFit::Detail::RooSimNLL; namespace { @@ -357,7 +358,7 @@ std::unique_ptr createSimultaneousNLL(RooSimultaneous const &simPdf, } // Time to sum the NLLs - auto nll = std::make_unique("mynll", "mynll", nllTerms); + auto nll = std::make_unique("mynll", "mynll", nllTerms, simCat, true); nll->addOwnedComponents(std::move(nllTerms)); return nll; } diff --git a/roofit/roofitcore/src/RooNLLVarNew.cxx b/roofit/roofitcore/src/RooNLLVarNew.cxx index f5edc2205fb17..120f0ad69024a 100644 --- a/roofit/roofitcore/src/RooNLLVarNew.cxx +++ b/roofit/roofitcore/src/RooNLLVarNew.cxx @@ -25,14 +25,15 @@ computation times. #include "RooFit/Detail/RooNLLVarNew.h" -#include +#include #include +#include #include +#include +#include #include -#include #include #include -#include #include "RooFitImplHelpers.h" @@ -47,6 +48,7 @@ computation times. #include ClassImp(RooFit::Detail::RooNLLVarNew); +ClassImp(RooFit::Detail::RooSimNLL); namespace RooFit { namespace Detail { @@ -336,6 +338,57 @@ void RooNLLVarNew::finalizeResult(RooFit::EvalContext &ctx, ROOT::Math::KahanSum ctx.setOutputWithOffset(this, result, _offset); } +RooSimNLL::RooSimNLL(const char *name, const char *title, const RooArgSet &terms, RooAbsCategoryLValue const &indexCat, + bool channelMasking) + : RooAbsReal(name, title), _set("!set", "set of components", this), _mask("!mask", "set of masks", this) +{ + _set.addTyped(terms); + + if (channelMasking) { + for (auto const &catState : indexCat) { + std::string const &catName = catState.first; + std::string maskName = "mask_" + catName; + _mask.addOwned(std::make_unique(maskName.c_str(), maskName.c_str(), 0.0)); + } + } +} + +RooSimNLL::RooSimNLL(const RooSimNLL &other, const char *name) + : RooAbsReal(other, name), _set("!set", this, other._set), _mask("!mask", this, other._set) +{ +} + +double RooSimNLL::evaluate() const +{ + double sum(0); + const RooArgSet *nset = _set.nset(); + + std::size_t i = 0; + for (auto *comp : static_range_cast(_set)) { + if (_mask.empty() || static_cast(_mask[i])->getVal() == 0.0) { + sum += comp->getVal(nset); + } + ++i; + } + return sum; +} + +void RooSimNLL::doEval(RooFit::EvalContext &ctx) const +{ + double result = 0.; + for (std::size_t i = 0; i < _set.size(); ++i) { + if (_mask.empty() || ctx.at(_mask[i])[0] == 0.0) { + result += ctx.at(_set[i])[0]; + } + } + ctx.output()[0] = result; +} + +double RooSimNLL::defaultErrorLevel() const +{ + return static_cast(_set[0])->defaultErrorLevel(); +} + } // namespace Detail } // namespace RooFit