From da8485368b87cbcbbe9d5ea1727ec498e2974eb8 Mon Sep 17 00:00:00 2001 From: will-cern Date: Fri, 19 Jun 2026 18:39:45 +0100 Subject: [PATCH 1/5] [RF] Update xRooFit Update xRooFit again, to accommodate more features/improvements made in past week while working on the Higgs discovery workspaces. (cherry picked from commit d536fdf8d3ad857e26d1647638c7cdf5b7d13332) --- roofit/xroofit/CMakeLists.txt | 1 + .../xroofit/inc/RooFit/xRooFit/xRooNLLVar.h | 20 +- roofit/xroofit/inc/RooFit/xRooFit/xRooNode.h | 35 +- roofit/xroofit/src/PythonInterface.cxx | 30 + roofit/xroofit/src/PythonInterface.h | 18 + roofit/xroofit/src/xRooFit.cxx | 241 +++++--- roofit/xroofit/src/xRooFitVersion.h | 4 +- roofit/xroofit/src/xRooHypoSpace.cxx | 87 ++- roofit/xroofit/src/xRooNLLVar.cxx | 127 +++-- roofit/xroofit/src/xRooNode.cxx | 519 +++++++++++------- 10 files changed, 728 insertions(+), 354 deletions(-) create mode 100644 roofit/xroofit/src/PythonInterface.cxx create mode 100644 roofit/xroofit/src/PythonInterface.h diff --git a/roofit/xroofit/CMakeLists.txt b/roofit/xroofit/CMakeLists.txt index c7c0a2638ffd5..edf683fd1cd61 100644 --- a/roofit/xroofit/CMakeLists.txt +++ b/roofit/xroofit/CMakeLists.txt @@ -15,6 +15,7 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitXRooFit src/xRooNLLVar.cxx src/xRooNode.cxx src/xRooNode_interactive.cxx + src/PythonInterface.cxx DICTIONARY_OPTIONS "-writeEmptyRootPCM" DEPENDENCIES diff --git a/roofit/xroofit/inc/RooFit/xRooFit/xRooNLLVar.h b/roofit/xroofit/inc/RooFit/xRooFit/xRooNLLVar.h index 9887bb4797231..922518e055144 100644 --- a/roofit/xroofit/inc/RooFit/xRooFit/xRooNLLVar.h +++ b/roofit/xroofit/inc/RooFit/xRooFit/xRooNLLVar.h @@ -63,6 +63,7 @@ class xRooNLLVar : public std::shared_ptr { xValueWithError(const std::pair &in = {0, 0}) : std::pair(in) {} double value() const { return std::pair::first; } double error() const { return std::pair::second; } + std::string __repr__() const { return Form("%g +/- %g", value(), error()); } }; void Print(Option_t *opt = ""); @@ -405,8 +406,8 @@ class xRooNLLVar : public std::shared_ptr { double alt_value = std::numeric_limits::quiet_NaN(), const xRooFit::Asymptotics::PLLType &pllType = xRooFit::Asymptotics::Unknown); xRooHypoSpace hypoSpace(const char *parName, xRooFit::TestStatistic::Type tsType, int nPoints = 0, - double low = -std::numeric_limits::infinity(), - double high = std::numeric_limits::infinity(), + double low = std::numeric_limits::quiet_NaN(), + double high = std::numeric_limits::quiet_NaN(), double alt_value = std::numeric_limits::quiet_NaN()) { return hypoSpace(parName, nPoints, low, high, alt_value, xRooFit::Asymptotics::Unknown, tsType); @@ -511,11 +512,18 @@ class xRooNLLVar : public std::shared_ptr { bool kReuseNLL = true; }; +END_XROOFIT_NAMESPACE + +#ifndef XROOFIT_NAMESPACE_NAME +#ifdef XROOFIT_NAMESPACE +#define XROOFIT_NAMESPACE_NAME XROOFIT_NAMESPACE +#else +#define XROOFIT_NAMESPACE_NAME +#endif +#endif namespace cling { -std::string printValue(const xRooNLLVar::xValueWithError *val); -std::string printValue(const std::map *m); +std::string printValue(const XROOFIT_NAMESPACE_NAME::xRooNLLVar::xValueWithError *val); +std::string printValue(const std::map *m); } // namespace cling -END_XROOFIT_NAMESPACE - #endif // include guard diff --git a/roofit/xroofit/inc/RooFit/xRooFit/xRooNode.h b/roofit/xroofit/inc/RooFit/xRooFit/xRooNode.h index df9f877b2493e..ee7c213968ca1 100644 --- a/roofit/xroofit/inc/RooFit/xRooFit/xRooNode.h +++ b/roofit/xroofit/inc/RooFit/xRooFit/xRooNode.h @@ -197,7 +197,23 @@ class xRooNode : public TNamed, public std::vector> { return aa == bb; }; }; - auto begin() const -> xRooNodeIterator { return xRooNodeIterator(std::vector>::begin()); } + auto begin() const -> xRooNodeIterator + { + // this update is to resolve need for calling browse() before iterating, e.g. + // for(auto a : xRooNode(b).browse()) ... + // can now become the more natural: + // for(auto a : xRooNode(b)) ... + // but would still need e.g. xRooNode(s).browse().size() rather than xRooNode(s).size() which would give 0 for + // unpopulated case + static bool browseLock = + false; // need for blocking recursive calls to browse (since browse method uses the iterators) + if (!browseLock && get() && empty()) { + browseLock = true; + const_cast(*this).browse(); + browseLock = false; + } + return xRooNodeIterator(std::vector>::begin()); + } auto end() const -> xRooNodeIterator { return xRooNodeIterator(std::vector>::end()); } // needed in pyROOT to avoid it creating iterators that follow the 'get' to death @@ -312,6 +328,9 @@ class xRooNode : public TNamed, public std::vector> { xRooNode datasets() const; // datasets corresponding to this pdf (parent nodes that do observable selections automatically applied) + xRooNode parents() const; // the clients of this node + xRooNode args() const; // all the nodes underneath this node + xRooNode Replace(const xRooNode &node); // use to replace a node in the tree at the location of this node xRooNode Remove(const xRooNode &child); xRooNode @@ -326,6 +345,7 @@ class xRooNode : public TNamed, public std::vector> { xRooNode reduced(const std::string &range = "", bool invert = false) const; // return a node representing reduced version of this node, will use the SetRange to reduce if blank + xRooNode reduced(const std::function selector) const; // following versions are for the menu in the GUI /** @private */ @@ -537,10 +557,17 @@ class xRooNode : public TNamed, public std::vector> { ClassDefOverride(xRooNode, 0) }; +END_XROOFIT_NAMESPACE + +#ifndef XROOFIT_NAMESPACE_NAME +#ifdef XROOFIT_NAMESPACE +#define XROOFIT_NAMESPACE_NAME XROOFIT_NAMESPACE +#else +#define XROOFIT_NAMESPACE_NAME +#endif +#endif namespace cling { -std::string printValue(const xRooNode *val); +std::string printValue(const XROOFIT_NAMESPACE_NAME::xRooNode *val); } -END_XROOFIT_NAMESPACE - #endif // include guard diff --git a/roofit/xroofit/src/PythonInterface.cxx b/roofit/xroofit/src/PythonInterface.cxx new file mode 100644 index 0000000000000..f6c32462d9666 --- /dev/null +++ b/roofit/xroofit/src/PythonInterface.cxx @@ -0,0 +1,30 @@ +#include "./PythonInterface.h" + +#include + +namespace xPython { + + // https://docs.python.org/3.11/c-api/init.html#c.Py_IsInitialized + bool isPythonInitialized() { + using fn_t = int (*)(); + static auto f = reinterpret_cast(gSystem->DynFindSymbol("*", "Py_IsInitialized")); + return f && f(); + } + + // https://docs.python.org/3/c-api/sys.html#c.PySys_WriteStdout + void writeStdoutLine(const char *msg) { + using fn_t = void (*)(const char *, ...); + static auto f = reinterpret_cast(gSystem->DynFindSymbol("*", "PySys_WriteStdout")); + if (f) + f("%s\n", msg); + } + + // https://docs.python.org/3/c-api/sys.html#c.PySys_WriteStderr + void writeStderrLine(const char *msg) { + using fn_t = void (*)(const char *, ...); + static auto f = reinterpret_cast(gSystem->DynFindSymbol("*", "PySys_WriteStderr")); + if (f) + f("%s\n", msg); + } + +} // namespace xPython \ No newline at end of file diff --git a/roofit/xroofit/src/PythonInterface.h b/roofit/xroofit/src/PythonInterface.h new file mode 100644 index 0000000000000..d9585c181116a --- /dev/null +++ b/roofit/xroofit/src/PythonInterface.h @@ -0,0 +1,18 @@ +#ifndef xRooFit_PythonInterface_h +#define xRooFit_PythonInterface_h + +// Helper functions to use the Python C API without introducint a link-time +// dependency on libpython. These are intended to be used only when xRooFit is +// used from Python, where libpython is loaded by default and the required +// symbols can be found with gSystem. + +namespace xPython { + + bool isPythonInitialized(); + + // Don't call these other functions if Python is not initialized! + void writeStdoutLine(const char *msg); + void writeStderrLine(const char *msg); +} // namespace xPython + +#endif \ No newline at end of file diff --git a/roofit/xroofit/src/xRooFit.cxx b/roofit/xroofit/src/xRooFit.cxx index 9288df5a17795..0707cb107c4e7 100644 --- a/roofit/xroofit/src/xRooFit.cxx +++ b/roofit/xroofit/src/xRooFit.cxx @@ -64,6 +64,10 @@ #include "TROOT.h" #include "TBrowser.h" +#include "TEnv.h" + +#include "./PythonInterface.h" + BEGIN_XROOFIT_NAMESPACE std::shared_ptr xRooFit::sDefaultNLLOptions = nullptr; @@ -222,9 +226,8 @@ xRooFit::generateFrom(RooAbsPdf &pdf, const RooFitResult &_fr, bool expected, in !(cClass && strcmp(cClass->GetName(), "SimpleGaussianConstraint") == 0)) { TString className = (cClass) ? cClass->GetName() : "undefined"; oocoutW((TObject *)nullptr, Generation) - << "xRooFit::generateFrom : constraint term " << thePdf->GetName() - << " of type " << className << " is a non-supported type - result might be not correct " - << std::endl; + << "xRooFit::generateFrom : constraint term " << thePdf->GetName() << " of type " + << className << " is a non-supported type - result might be not correct " << std::endl; } // in case of a Poisson constraint make sure the rounding is not set @@ -514,7 +517,9 @@ std::shared_ptr xRooFit::defaultFitConfig() extraOpts->SetValue("OptimizeConst", 2); // if 0 will disable constant term optimization and cache-and-track of the // NLL. 1 = just caching, 2 = cache and track #if ROOT_VERSION_CODE >= ROOT_VERSION(6, 29, 00) - extraOpts->SetValue("StrategySequence", "0s01s12s2s3m"); + extraOpts->SetValue( + "StrategySequence", + "0s01s12s2r2s23"); // 'm' would indicate use migradImproved from minuit v1. Dropped from default for 6.40 extraOpts->SetValue("HesseStrategySequence", "23"); #else extraOpts->SetValue("StrategySequence", "0s01s12s2m"); @@ -541,6 +546,29 @@ ROOT::Math::IOptions *xRooFit::defaultFitConfigOptions() return const_cast(defaultFitConfig()->MinimizerOptions().ExtraOptions()); } +void printCerr(const char *msg) +{ + if (xPython::isPythonInitialized()) { + xPython::writeStderrLine(msg); + // if (PyObject *sys_stdout = PySys_GetObject("stderr"); sys_stdout != nullptr) { + // Py_XDECREF(PyObject_CallMethod(sys_stdout, "flush", nullptr)); + // } + } else { + std::cerr << msg << std::endl; + } +} +void printCout(const char *msg) +{ + if (xPython::isPythonInitialized()) { + xPython::writeStdoutLine(msg); + // if (PyObject *sys_stdout = PySys_GetObject("stdout"); sys_stdout != nullptr) { + // Py_XDECREF(PyObject_CallMethod(sys_stdout, "flush", nullptr)); + // } + } else { + std::cout << msg << std::endl; + } +} + class ProgressMonitor : public RooAbsReal { public: void (*oldHandlerr)(int) = nullptr; @@ -549,7 +577,7 @@ class ProgressMonitor : public RooAbsReal { static void interruptHandler(int signum) { if (signum == SIGINT) { - std::cout << "Minimization interrupted ... will exit as soon as possible" << std::endl; + printCout("Minimization interrupted ... will exit as soon as possible"); // TODO: create a global mutex for this fInterrupt = true; } else { @@ -625,14 +653,14 @@ class ProgressMonitor : public RooAbsReal { s.Reset(); std::stringstream sout; - sout << (counter) << ") (" << evalRate << "Hz) " << TDatime().AsString(); + sout << TDatime().AsString() << ":(" << (counter) << "|" << evalRate << "Hz)"; if (!fState.empty()) sout << " : " << fState; if (counter2) { // doing a hesse step, estimate progress based on evaluations int nRequired = prevPars.size(); if (nRequired > 1) { - nRequired *= (nRequired-1); + nRequired *= (nRequired - 1); nRequired /= 2; // since only need to do the a 'triangle' of the hessian matrix if (fState == "Hesse3") { nRequired *= 4; @@ -640,7 +668,7 @@ class ProgressMonitor : public RooAbsReal { sout << " (~" << int(100.0 * (counter - counter2) / nRequired) << "%)"; } } - sout << " : " << minVal << " Delta = " << (minVal - prevMin); + sout << " : " << minVal << " Delta=" << (minVal - prevMin); if (minVal < prevMin) { sout << " : "; // compare minPars and prevPars, print biggest deltas @@ -685,7 +713,7 @@ class ProgressMonitor : public RooAbsReal { } gSystem->ProcessEvents(); } - std::cerr << sout.str() << std::endl; + printCerr(sout.str().c_str()); // std::cerr << sout.str() << std::endl; prevMin = minVal; prevCounter = counter; @@ -699,10 +727,11 @@ class ProgressMonitor : public RooAbsReal { mutable int counter = 0; int counter2 = 0; // used to estimate progress of a Hesse calculation -private: - RooRealProxy fFunc; mutable double minVal = std::numeric_limits::infinity(); mutable double prevMin = std::numeric_limits::infinity(); + +private: + RooRealProxy fFunc; mutable RooArgList minPars; mutable RooArgList prevPars; mutable int prevCounter = 0; @@ -948,13 +977,14 @@ std::shared_ptr xRooFit::minimize(RooAbsReal &nll, if (!out) { int strategy = fitConfig.MinimizerOptions().Strategy(); // Note: AsymptoticCalculator enforces not less than 1 on tolerance - should we do so too? - if (_progress) { + if (_progress && printLevel >= -2) { _nll = new ProgressMonitor(*_nll, _progress); ProgressMonitor::fInterrupt = false; } auto logger = (logSize > 0) ? std::make_unique(logs, logSize) : nullptr; - RooMinimizer _minimizer(*_nll); - _minimizer.fitter()->Config() = fitConfig; + std::unique_ptr _minimizerPtr = std::make_unique(*_nll); + RooMinimizer *_minimizer = _minimizerPtr.get(); + _minimizer->fitter()->Config() = fitConfig; // if(fitConfig.MinimizerOptions().ExtraOptions()) { // //for loading hesse options // double a; @@ -966,34 +996,35 @@ std::shared_ptr xRooFit::minimize(RooAbsReal &nll, // } // } - bool autoMaxCalls = (_minimizer.fitter()->Config().MinimizerOptions().MaxFunctionCalls() == 0); + bool autoMaxCalls = (_minimizer->fitter()->Config().MinimizerOptions().MaxFunctionCalls() == 0); if (autoMaxCalls) { - _minimizer.fitter()->Config().MinimizerOptions().SetMaxFunctionCalls( + _minimizer->fitter()->Config().MinimizerOptions().SetMaxFunctionCalls( 500 * floatPars->size() * floatPars->size()); // hesse requires O(N^2) function calls } - if (_minimizer.fitter()->Config().MinimizerOptions().MaxIterations() == 0) { - _minimizer.fitter()->Config().MinimizerOptions().SetMaxIterations(500 * floatPars->size()); + if (_minimizer->fitter()->Config().MinimizerOptions().MaxIterations() == 0) { + _minimizer->fitter()->Config().MinimizerOptions().SetMaxIterations(500 * floatPars->size()); } - bool hesse = _minimizer.fitter()->Config().ParabErrors(); - _minimizer.fitter()->Config().SetParabErrors( + std::unique_ptr floatInitVals(floatPars->snapshot()); + bool hesse = _minimizer->fitter()->Config().ParabErrors(); + _minimizer->fitter()->Config().SetParabErrors( false); // turn "off" so can run hesse as a separate step, appearing in status - _minimizer.fitter()->Config().SetMinosErrors(false); - _minimizer.fitter()->Config().SetUpdateAfterFit(true); // note: seems to always take effect + _minimizer->fitter()->Config().SetMinosErrors(false); + _minimizer->fitter()->Config().SetUpdateAfterFit(true); // note: seems to always take effect std::vector> statusHistory; // gCurrentSampler = this; // gOldHandlerr = signal(SIGINT,toyInterruptHandlerr); - TString actualFirstMinimizer = _minimizer.fitter()->Config().MinimizerType(); + TString actualFirstMinimizer = _minimizer->fitter()->Config().MinimizerType(); int status = 0; int constOptimize = 2; - _minimizer.fitter()->Config().MinimizerOptions().ExtraOptions()->GetValue("OptimizeConst", constOptimize); + _minimizer->fitter()->Config().MinimizerOptions().ExtraOptions()->GetValue("OptimizeConst", constOptimize); if (constOptimize) { - _minimizer.optimizeConst(constOptimize); + _minimizer->optimizeConst(constOptimize); // for safety force a refresh of the cache (and tracking) in the nll // DO NOT do a ConfigChange ... this is just a deactivate-reactivate of caching // but it seems like doing this breaks the const optimization and function is badly behaved @@ -1011,8 +1042,8 @@ std::shared_ptr xRooFit::minimize(RooAbsReal &nll, } int sIdx = -1; - TString minim = _minimizer.fitter()->Config().MinimizerType(); - TString algo = _minimizer.fitter()->Config().MinimizerAlgoType(); + TString minim = _minimizer->fitter()->Config().MinimizerType(); + TString algo = _minimizer->fitter()->Config().MinimizerAlgoType(); if (minim == "Minuit2") { if (strategy == -1) { sIdx = 0; @@ -1038,22 +1069,41 @@ std::shared_ptr xRooFit::minimize(RooAbsReal &nll, algo = "Scan"; } else if (m_strategy(sIdx) == 'h') { break; // jumping straight to a hesse evaluation + } else if (m_strategy(sIdx) == 'r') { + // reset minimizer + resetMinimization: + tries = 0; + *floatPars = *floatInitVals; // resets floats + std::unique_ptr _minimizerPtr2 = std::make_unique(*_nll); + auto initPars = _minimizerPtr2->fitter()->Config().ParamsSettings(); + auto initOpts = _minimizerPtr2->fitter()->Config().MinimizerOptions(); + _minimizerPtr2->fitter()->Config() = _minimizer->fitter()->Config(); + _minimizerPtr2->fitter()->Config().SetParamsSettings({}); // clears param settings + _minimizerPtr = std::move(_minimizerPtr2); + _minimizer = _minimizerPtr.get(); + sIdx++; + statusHistory.emplace_back("Reset", 0); + if (auto fff = dynamic_cast(_nll); fff) { + fff->counter2 = 0; // may have become non-zero if progressed to hesse and then resumed + fff->prevMin = fff->minVal; // reset minimum + } + continue; } else { strategy = int(m_strategy(sIdx) - '0'); - _minimizer.setStrategy(strategy); + _minimizer->setStrategy(strategy); minim = "Minuit2"; algo = "Migrad"; } if (auto fff = dynamic_cast(_nll); fff) { - fff->fState = minim + algo + std::to_string(_minimizer.fitter()->Config().MinimizerOptions().Strategy()); + fff->fState = minim + algo + std::to_string(_minimizer->fitter()->Config().MinimizerOptions().Strategy()); } try { - status = _minimizer.minimize(minim, algo); + status = _minimizer->minimize(minim, algo); } catch (const std::exception &e) { std::cerr << "Exception while minimizing: " << e.what() << std::endl; } - if (first && actualFirstMinimizer != _minimizer.fitter()->Config().MinimizerType()) - actualFirstMinimizer = _minimizer.fitter()->Config().MinimizerType(); + if (first && actualFirstMinimizer != _minimizer->fitter()->Config().MinimizerType()) + actualFirstMinimizer = _minimizer->fitter()->Config().MinimizerType(); first = false; tries++; @@ -1063,13 +1113,13 @@ std::shared_ptr xRooFit::minimize(RooAbsReal &nll, } // RooMinimizer loses the useful status code, so here we will override it - status = _minimizer.fitter() + status = _minimizer->fitter() ->Result() .Status(); // note: Minuit failure is status code 4, minuit2 that is edm above max - minim = _minimizer.fitter()->Config().MinimizerType(); // may have changed value - statusHistory.emplace_back(_minimizer.fitter()->Config().MinimizerType() + - _minimizer.fitter()->Config().MinimizerAlgoType() + - std::to_string(_minimizer.fitter()->Config().MinimizerOptions().Strategy()), + minim = _minimizer->fitter()->Config().MinimizerType(); // may have changed value + statusHistory.emplace_back(_minimizer->fitter()->Config().MinimizerType() + + _minimizer->fitter()->Config().MinimizerAlgoType() + + std::to_string(_minimizer->fitter()->Config().MinimizerOptions().Strategy()), status); if (status % 1000 == 0) break; // fit was good @@ -1077,15 +1127,15 @@ std::shared_ptr xRooFit::minimize(RooAbsReal &nll, if (status == 4 && minim != "Minuit") { if (printLevel >= -1) { Warning("fitTo", "%s Hit max function calls of %d", fitName.Data(), - _minimizer.fitter()->Config().MinimizerOptions().MaxFunctionCalls()); + _minimizer->fitter()->Config().MinimizerOptions().MaxFunctionCalls()); } if (autoMaxCalls) { if (printLevel >= -1) Warning("fitTo", "will try doubling this"); - _minimizer.fitter()->Config().MinimizerOptions().SetMaxFunctionCalls( - _minimizer.fitter()->Config().MinimizerOptions().MaxFunctionCalls() * 2); - _minimizer.fitter()->Config().MinimizerOptions().SetMaxIterations( - _minimizer.fitter()->Config().MinimizerOptions().MaxIterations() * 2); + _minimizer->fitter()->Config().MinimizerOptions().SetMaxFunctionCalls( + _minimizer->fitter()->Config().MinimizerOptions().MaxFunctionCalls() * 2); + _minimizer->fitter()->Config().MinimizerOptions().SetMaxIterations( + _minimizer->fitter()->Config().MinimizerOptions().MaxIterations() * 2); continue; } } @@ -1094,11 +1144,13 @@ std::shared_ptr xRooFit::minimize(RooAbsReal &nll, // specified Also note that if fits are failing because of edm over max, it can be a good idea to activate the // Offset option when building nll if (printLevel >= -1) { - Warning("fitTo", "%s %s%s Status=%d (edm=%f, tol=%f, strat=%d), tries=#%d...", fitName.Data(), - _minimizer.fitter()->Config().MinimizerType().c_str(), - _minimizer.fitter()->Config().MinimizerAlgoType().c_str(), status, - _minimizer.fitter()->Result().Edm(), _minimizer.fitter()->Config().MinimizerOptions().Tolerance(), - _minimizer.fitter()->Config().MinimizerOptions().Strategy(), tries); + printCerr(TString::Format("Warning: %s %s%s Status=%d (edm=%f, tol=%f, strat=%d), tries=#%d...", + fitName.Data(), _minimizer->fitter()->Config().MinimizerType().c_str(), + _minimizer->fitter()->Config().MinimizerAlgoType().c_str(), status, + _minimizer->fitter()->Result().Edm(), + _minimizer->fitter()->Config().MinimizerOptions().Tolerance(), + _minimizer->fitter()->Config().MinimizerOptions().Strategy(), tries) + .Data()); } // decide what to do next based on strategy sequence @@ -1109,6 +1161,7 @@ std::shared_ptr xRooFit::minimize(RooAbsReal &nll, tries--; sIdx++; } + auto mini_sIdx = sIdx; /* Minuit2 status codes: * status = 0 : OK @@ -1131,16 +1184,16 @@ std::shared_ptr xRooFit::minimize(RooAbsReal &nll, // Note that if strategy>=2 or (strategy=1 and Dcovar>0.05) then hesse will be forced to be run (see // VariadicMetricBuilder) So only in Strategy=0 can you skip hesse (even if SetParabErrors false). - int miniStrat = _minimizer.fitter()->Config().MinimizerOptions().Strategy(); + int miniStrat = _minimizer->fitter()->Config().MinimizerOptions().Strategy(); double dCovar = std::numeric_limits::quiet_NaN(); - // if(auto _minuit2 = dynamic_cast(_minimizer.fitter()->GetMinimizer()); + // if(auto _minuit2 = dynamic_cast(_minimizer->fitter()->GetMinimizer()); // _minuit2 && _minuit2->fMinimum) { // dCovar = _minuit2->fMinimum->Error().Dcovar(); // } // only do hesse if was a valid min and not full accurate cov matrix already (can happen if e.g. ran strat2) if (hesse && m_hessestrategy.Length() != 0 && - (m_strategy(sIdx) == 'h' || (_minimizer.fitter()->Result().IsValid()))) { + (m_strategy(sIdx) == 'h' || (_minimizer->fitter()->Result().IsValid()))) { // Note: minima where the covariance was made posdef are deemed 'valid' ... @@ -1148,8 +1201,8 @@ std::shared_ptr xRooFit::minimize(RooAbsReal &nll, // transformed interesting note: error on pars before hesse can be significantly smaller than after hesse ... // what is the pre-hesse error corresponding to? - corresponds to approximation of covariance matrix calculated // with iterative method - /*auto parSettings = _minimizer.fitter()->Config().ParamsSettings(); - for (auto &ss : _minimizer.fitter()->Config().ParamsSettings()) { + /*auto parSettings = _minimizer->fitter()->Config().ParamsSettings(); + for (auto &ss : _minimizer->fitter()->Config().ParamsSettings()) { ss.RemoveLimits(); } @@ -1159,8 +1212,8 @@ std::shared_ptr xRooFit::minimize(RooAbsReal &nll, v->removeRange(); }*/ - // std::cout << "nIterations = " << _minimizer.fitter()->GetMinimizer()->NIterations() << std::endl; - // std::cout << "covQual before hesse = " << _minimizer.fitter()->GetMinimizer()->CovMatrixStatus() << + // std::cout << "nIterations = " << _minimizer->fitter()->GetMinimizer()->NIterations() << std::endl; + // std::cout << "covQual before hesse = " << _minimizer->fitter()->GetMinimizer()->CovMatrixStatus() << // std::endl; sIdx = -1; if (hesseStrategy == -1) { @@ -1179,7 +1232,7 @@ std::shared_ptr xRooFit::minimize(RooAbsReal &nll, if (strategy == 2 && hesseStrategy == 2) { // don't repeat hesse if strategy=2 and hesseStrategy=2, and the matrix was valid - if (_minimizer.fitter()->GetMinimizer()->CovMatrixStatus() == 3) { + if (_minimizer->fitter()->GetMinimizer()->CovMatrixStatus() == 3) { break; } if (sIdx >= m_hessestrategy.Length() - 1) { @@ -1189,22 +1242,24 @@ std::shared_ptr xRooFit::minimize(RooAbsReal &nll, continue; } - _minimizer.fitter()->Config().MinimizerOptions().SetStrategy(hesseStrategy); - // const_cast(_minimizer.fitter()->Config().MinimizerOptions().ExtraOptions())->SetValue("HessianStepTolerance",0.1); - // const_cast(_minimizer.fitter()->Config().MinimizerOptions().ExtraOptions())->SetValue("HessianG2Tolerance",0.02); + _minimizer->fitter()->Config().MinimizerOptions().SetStrategy(hesseStrategy); + // const_cast(_minimizer->fitter()->Config().MinimizerOptions().ExtraOptions())->SetValue("HessianStepTolerance",0.1); + // const_cast(_minimizer->fitter()->Config().MinimizerOptions().ExtraOptions())->SetValue("HessianG2Tolerance",0.02); if (auto fff = dynamic_cast(_nll); fff) { - fff->fState = TString::Format("Hesse%d", _minimizer.fitter()->Config().MinimizerOptions().Strategy()); + fff->fState = TString::Format("Hesse%d", _minimizer->fitter()->Config().MinimizerOptions().Strategy()); fff->counter2 = fff->counter; + fff->prevMin = + fff->minVal; // reset minimum when change to hesse. Helps see if hesse eval gives new lower values } //_nll->getVal(); // for reasons I dont understand, if nll evaluated before hesse call the edm is smaller? - // and also becomes WRONG :-S - // auto _status = (_minimizer.fitter()->CalculateHessErrors()) ? _minimizer.fitter()->Result().Status() : + // auto _status = (_minimizer->fitter()->CalculateHessErrors()) ? _minimizer->fitter()->Result().Status() : // -1; - auto _status = _minimizer.hesse(); // note: I have seen that you can get 'full covariance quality' without - // running hesse ... is that expected? + auto _status = _minimizer->hesse(); // note: I have seen that you can get 'full covariance quality' without + // running hesse ... is that expected? // note: hesse status will be -1 if hesse failed (no covariance matrix) // otherwise the status appears to be whatever was the status before // note that hesse succeeds even if the cov matrix it calculates is forced pos def. Failure is only @@ -1219,39 +1274,48 @@ std::shared_ptr xRooFit::minimize(RooAbsReal &nll, v->removeRange("backup"); } } - _minimizer.fitter()->Config().SetParamsSettings(parSettings);*/ + _minimizer->fitter()->Config().SetParamsSettings(parSettings);*/ - /*for (auto &ss : _minimizer.fitter()->Config().ParamsSettings()) { + /*for (auto &ss : _minimizer->fitter()->Config().ParamsSettings()) { if( ss.HasLowerLimit() || ss.HasUpperLimit() ) std::cout << ss.Name() << " limit restored " << ss.LowerLimit() << " - " << ss.UpperLimit() << std::endl; }*/ statusHistory.push_back(std::pair( - TString::Format("Hesse%d", _minimizer.fitter()->Config().MinimizerOptions().Strategy()), _status)); + TString::Format("Hesse%d", _minimizer->fitter()->Config().MinimizerOptions().Strategy()), _status)); if (auto fff = dynamic_cast(_nll); fff && fff->fInterrupt) { delete _nll; throw std::runtime_error("Keyboard interrupt while hesse calculating"); } - if ((_status != 0 || _minimizer.fitter()->GetMinimizer()->CovMatrixStatus() != 3) && status == 0 && + if ((_status != 0 || _minimizer->fitter()->GetMinimizer()->CovMatrixStatus() != 3) && status == 0 && printLevel >= -1) { - Warning("fitTo", "%s hesse status is %d, covQual=%d", fitName.Data(), _status, - _minimizer.fitter()->GetMinimizer()->CovMatrixStatus()); - } - - if (sIdx >= m_hessestrategy.Length() - 1) { - break; // run out of strategies to try, stop + printCerr(TString::Format("Warning: %s hesse status is %d, covQual=%d", fitName.Data(), _status, + _minimizer->fitter()->GetMinimizer()->CovMatrixStatus()) + .Data()); } - if (_status == 0 && _minimizer.fitter()->GetMinimizer()->CovMatrixStatus() == 3) { + if (_status == 0 && _minimizer->fitter()->GetMinimizer()->CovMatrixStatus() == 3) { // covariance is valid! break; } else if (_status == 0) { // set the statusHistory to the cov status, since that's more informative - statusHistory.back().second = _minimizer.fitter()->GetMinimizer()->CovMatrixStatus(); + statusHistory.back().second = _minimizer->fitter()->GetMinimizer()->CovMatrixStatus(); + } + + if (sIdx >= m_hessestrategy.Length() - 1) { + break; // run out of strategies to try, stop } + sIdx++; } // end of hesse attempt loop + // experimental feature to resume fits invalidated by hesse + if (gEnv->GetValue("XRooFit.ResumeInvalidFits", false) && + (statusHistory.back().second == 1 || statusHistory.back().second == 2) && + mini_sIdx < (m_strategy.Length() - 1)) { + sIdx = mini_sIdx; + goto resetMinimization; + } } // call minos if requested on any parameters @@ -1260,18 +1324,27 @@ std::shared_ptr xRooFit::minimize(RooAbsReal &nll, if (auto fff = dynamic_cast(_nll); fff) { fff->fState = "Minos"; fff->counter2 = 0; + fff->prevMin = fff->minVal; } - auto _status = _minimizer.minos(*mpars); + auto _status = _minimizer->minos(*mpars); statusHistory.push_back(std::pair("Minos", _status)); } } // DO NOT DO THIS - seems to mess with the NLL function in a way that breaks the cache - reactivating wont fix - // if(constOptimize) { _minimizer.optimizeConst(0); } // doing this because saw happens in RooAbsPdf::minimizeNLL + // if(constOptimize) { _minimizer->optimizeConst(0); } // doing this because saw happens in RooAbsPdf::minimizeNLL // method // signal(SIGINT,gOldHandlerr); - out = std::unique_ptr{_minimizer.save(fitName, resultTitle)}; + out = std::unique_ptr{_minimizer->save(fitName, resultTitle)}; + + // if the final result is not valid, RooFit will mark the status as -1. But we should override that with the + // more informative status of the last fit ... + // for safety, will only override if status of last result is not 0 (don't want to report success when actually + // bad) + if (out->status() == -1 && !_minimizer->fitter()->Result().IsValid() && _minimizer->fitter()->Result().Status()) { + out->setStatus(_minimizer->fitter()->Result().Status()); + } // if status is 0 (min succeeded) but the covQual isn't fully accurate but requested hesse, reflect that in the // status @@ -1284,13 +1357,15 @@ std::shared_ptr xRooFit::minimize(RooAbsReal &nll, } } - if (miniStrat < _minimizer.fitter()->Config().MinimizerOptions().Strategy() && hesse && - out->edm() > _minimizer.fitter()->Config().MinimizerOptions().Tolerance() * 1e-3 && out->status() != 3) { + if (printLevel >= -2 && miniStrat < _minimizer->fitter()->Config().MinimizerOptions().Strategy() && hesse && + out->edm() > _minimizer->fitter()->Config().MinimizerOptions().Tolerance() * 1e-2 && out->status() != 3) { // hesse may have updated edm by using a better strategy than used in the minimization // so print a warning about this - std::cerr << "Warning: post-Hesse edm " << out->edm() << " greater than allowed by tolerance " - << _minimizer.fitter()->Config().MinimizerOptions().Tolerance() * 1e-3 - << ", consider increasing minimization strategy" << std::endl; + std::stringstream ss; + ss << "Warning: post-Hesse edm " << out->edm() + << " > 10xMaxEDM (MaxEDM=" << _minimizer->fitter()->Config().MinimizerOptions().Tolerance() * 1e-3 + << "). Consider increasing your minimization strategy"; + printCerr(ss.str().c_str()); // Dec24: As this is a new warning, will not update status code for now, so edm will be large // but in the future we should probably update the code to 3 so that users don't miss this warning. // out->setStatus(3); // edm above max @@ -1399,7 +1474,7 @@ std::shared_ptr xRooFit::minimize(RooAbsReal &nll, fitConfig.MinimizerOptions().SetMinimizerType(actualFirstMinimizer); } - if (_progress) { + if (_progress && printLevel >= -2) { delete _nll; } } diff --git a/roofit/xroofit/src/xRooFitVersion.h b/roofit/xroofit/src/xRooFitVersion.h index cf46130bb4d64..27854f6e95f07 100644 --- a/roofit/xroofit/src/xRooFitVersion.h +++ b/roofit/xroofit/src/xRooFitVersion.h @@ -12,5 +12,5 @@ #pragma once -#define GIT_COMMIT_HASH "v0.0.2" -#define GIT_COMMIT_DATE "2026-04-28 15:01:00 +0100" +#define GIT_COMMIT_HASH "v0.0.4" +#define GIT_COMMIT_DATE "2026-06-19 16:32:00 +0100" diff --git a/roofit/xroofit/src/xRooHypoSpace.cxx b/roofit/xroofit/src/xRooHypoSpace.cxx index f016ee2ac791f..e8a0714f4f898 100644 --- a/roofit/xroofit/src/xRooHypoSpace.cxx +++ b/roofit/xroofit/src/xRooHypoSpace.cxx @@ -37,6 +37,8 @@ #include "RooStats/HypoTestInverterResult.h" #include "TEnv.h" +#include "./PythonInterface.h" + BEGIN_XROOFIT_NAMESPACE xRooNLLVar::xRooHypoSpace::xRooHypoSpace(const char *name, const char *title) @@ -316,16 +318,14 @@ int xRooNLLVar::xRooHypoSpace::scan(const char *type, size_t nPoints, double low } if (/*high < low ||*/ (high == low && nPoints != 1)) { - // take from parameter + // take from parameter (will be either the defaults from the construction of the hypoSpace or whatever last scan + // was low = p->getMin("scan"); high = p->getMax("scan"); } if (!std::isnan(low) && !std::isnan(high) && !(std::isinf(low) && std::isinf(high))) { p->setRange("scan", std::min(low, high), std::max(low, high)); } - if (p->hasRange("scan")) { - ::Info("xRooHypoSpace::scan", "Using %s scan range: %g - %g", p->GetName(), p->getMin("scan"), p->getMax("scan")); - } bool doObs = false; for (auto nSigma : nSigmas) { @@ -403,6 +403,10 @@ int xRooNLLVar::xRooHypoSpace::scan(const char *type, size_t nPoints, double low if (nPoints == 0) { // automatic scan if (sType.Contains("cls")) { + if (p->hasRange("scan")) { + ::Info("xRooHypoSpace::scan", "Using %s scan range: %g - %g", p->GetName(), p->getMin("scan"), + p->getMax("scan")); + } for (double nSigma : nSigmas) { xValueWithError res(std::make_pair(0., 0.)); if (std::isnan(nSigma)) { @@ -428,8 +432,8 @@ int xRooNLLVar::xRooHypoSpace::scan(const char *type, size_t nPoints, double low out += 1; } else { // fit was ok, so use the values to determine an appropriate range - low = std::max(low, back().mu_hat().getVal()-back().mu_hat().getError()*3); - high = std::min(high, back().mu_hat().getVal()+back().mu_hat().getError()*3); + low = std::max(low, back().mu_hat().getVal() - back().mu_hat().getError() * 3); + high = std::min(high, back().mu_hat().getVal() + back().mu_hat().getError() * 3); nPoints = 20; double step = (high - low) / (nPoints - 1); for (size_t i = 0; i < nPoints; i++) { @@ -438,7 +442,6 @@ int xRooNLLVar::xRooHypoSpace::scan(const char *type, size_t nPoints, double low if (back().status() != 0) out += 1; } - } } else { throw std::runtime_error(TString::Format("Automatic scanning not yet supported for %s", type)); @@ -642,7 +645,15 @@ xRooNLLVar::xRooHypoPoint &xRooNLLVar::xRooHypoSpace::AddPoint(const char *coord } coordString.erase(coordString.end() - 1); - ::Info("xRooHypoSpace::AddPoint", "Added new point @ %s", coordString.c_str()); + if (xPython::isPythonInitialized()) { + auto s = TString::Format("Info in : Added new point @ %s", coordString.c_str()); + xPython::writeStdoutLine(s.Data()); + // if (PyObject *sys_stdout = PySys_GetObject("stdout"); sys_stdout != nullptr) { + // Py_XDECREF(PyObject_CallMethod(sys_stdout, "flush", nullptr)); + // } + } else { + ::Info("xRooHypoSpace::AddPoint", "Added new point @ %s", coordString.c_str()); + } return emplace_back(out); } @@ -899,8 +910,8 @@ std::shared_ptr xRooNLLVar::xRooHypoSpace::graph( out->GetListOfFunctions()->Add(x, "F"); x = out->Clone("down"); x->SetBit(kCanDelete); - // dynamic_cast(x)->SetFillColor((nSigma==2) ? kYellow : kGreen); - // dynamic_cast(x)->SetFillStyle(1001); + dynamic_cast(x)->SetFillColor(kBlack); + dynamic_cast(x)->SetFillStyle(nSigma == 2 ? 3005 : 3004); out->GetListOfFunctions()->Add(x, "F"); } if (sOpt.Contains("ts")) { @@ -1006,12 +1017,22 @@ std::shared_ptr xRooNLLVar::xRooHypoSpace::graph( for (int i = 0; i < out->GetN(); i++) { if (i < out->GetN() - nPointsDown) { up->SetPoint(up->GetN(), out->GetPointX(i), out->GetPointY(i) + out->GetErrorY(i) * (above ? 1. : -1.)); - down->SetPoint(down->GetN(), out->GetPointX(i), out->GetPointY(i) - out->GetErrorY(i) * (above ? 1. : -1.)); + //down->SetPoint(down->GetN(), out->GetPointX(i), out->GetPointY(i) - out->GetErrorY(i) * (above ? 1. : -1.)); } else { - up->SetPoint(up->GetN(), out->GetPointX(i), out->GetPointY(i) - out->GetErrorY(i) * (above ? 1. : -1.)); + //up->SetPoint(up->GetN(), out->GetPointX(i), out->GetPointY(i) - out->GetErrorY(i) * (above ? 1. : -1.)); down->SetPoint(down->GetN(), out->GetPointX(i), out->GetPointY(i) + out->GetErrorY(i) * (above ? 1. : -1.)); } } + // now go back round in reverse + for (int i = out->GetN()-1; i >= 0; i--) { + if (i < out->GetN() - nPointsDown) { + up->SetPoint(up->GetN(), out->GetPointX(i), out->GetPointY(i) - out->GetErrorY(i) * (above ? 1. : -1.)); + //down->SetPoint(down->GetN(), out->GetPointX(i), out->GetPointY(i) - out->GetErrorY(i) * (above ? 1. : -1.)); + } else { + //up->SetPoint(up->GetN(), out->GetPointX(i), out->GetPointY(i) - out->GetErrorY(i) * (above ? 1. : -1.)); + down->SetPoint(down->GetN(), out->GetPointX(i), out->GetPointY(i) - out->GetErrorY(i) * (above ? 1. : -1.)); + } + } } if (visualize) { @@ -1099,8 +1120,8 @@ std::shared_ptr xRooNLLVar::xRooHypoSpace::graphs(const char *opt) // out->GetListOfFunctions()->Add(out->GetHistogram()->Clone(".axis"),"sameaxis"); // redraw axis for (auto g : *out->GetListOfGraphs()) { - if (auto o = dynamic_cast(g)->GetListOfFunctions()->FindObject("down")) { - leg->AddEntry(o, "", "F"); + if (dynamic_cast(g)->GetListOfFunctions()->FindObject("down")) { + leg->AddEntry(g, "", "F"); } else { leg->AddEntry(g, "", "LPE"); } @@ -1154,7 +1175,8 @@ std::shared_ptr xRooNLLVar::xRooHypoSpace::graphs(const char *opt) auto gra2 = static_cast(out->DrawClone("A")); gra2->SetBit(kCanDelete); if (sOpt.Contains("pcls") || sOpt.Contains("pnull")) { - gra2->GetHistogram()->SetMinimum(1e-6); + gra2->SetMinimum(1e-6); + gra2->SetMaximum(1); } if (gPad) { gPad->RedrawAxis(); @@ -1268,11 +1290,13 @@ xRooNLLVar::xRooHypoSpace::findlimit(const char *opt, double relUncert, unsigned if (!gPad) gra->Draw(); // in 6.28 DrawClone wont make the gPad defined :( ... so Draw then clear and Draw Clone gPad->Clear(); - gra->DrawClone("A")->SetBit(kCanDelete); + auto gra2 = static_cast(gra->DrawClone("A")); + gra2->SetBit(kCanDelete); + gra2->SetMinimum(1e-9); + gra2->SetMaximum(1); gPad->RedrawAxis(); - gra->GetHistogram()->SetMinimum(1e-9); - gra->GetHistogram()->GetYaxis()->SetRangeUser(1e-9, 1); - gPad->Modified(); + gPad->GetCanvas()->Paint(); + gPad->GetCanvas()->Update(); #if ROOT_VERSION_CODE >= ROOT_VERSION(6, 30, 00) gPad->GetCanvas()->ResetUpdated(); // stops previous canvas being replaced in a jupyter notebook #endif @@ -1298,7 +1322,7 @@ xRooNLLVar::xRooHypoSpace::findlimit(const char *opt, double relUncert, unsigned if (!gr || gr->GetN() < 1) { if (maxTries == 0 || std::isnan(AddPoint(TString::Format("%s=%g", v->GetName(), muMin)).getVal(sOpt).first)) { // first point failed ... give up - ::Error("findlimit", "Problem evaluating %s @ %s=%g", sOpt.Data(), v->GetName(), muMin); + ::Error("findlimit", "Problem evaluating First Point %s @ %s=%g", sOpt.Data(), v->GetName(), muMin); return std::pair(std::numeric_limits::quiet_NaN(), 0.); } gr.reset(); @@ -1334,7 +1358,8 @@ xRooNLLVar::xRooHypoSpace::findlimit(const char *opt, double relUncert, unsigned if (maxTries == 0 || std::isnan(AddPoint(TString::Format("%s=%g", v->GetName(), nextPoint)).getVal(sOpt).first)) { // second point failed ... give up - ::Error("xRooHypoSpace::findlimit", "Problem evaluating %s @ %s=%g", sOpt.Data(), v->GetName(), nextPoint); + ::Error("xRooHypoSpace::findlimit", "Problem evaluating Second Point %s @ %s=%g", sOpt.Data(), v->GetName(), + nextPoint); return std::pair(std::numeric_limits::quiet_NaN(), 0.); } gr.reset(); @@ -1396,14 +1421,15 @@ xRooNLLVar::xRooHypoSpace::findlimit(const char *opt, double relUncert, unsigned // got here need a new point .... evaluate the estimated lim location +/- the relUncert (signed error takes care of // direction) + if (maxTries == 0) { + ::Warning("xRooHypoSpace::findlimit", "Reached max number of point evaluations"); + return lim; + } + ::Info("xRooHypoSpace::findlimit", "%s -- Testing new point @ %s=%g (delta=%g)", sOpt.Data(), v->GetName(), nextPoint, lim.second); - if (maxTries == 0 || std::isnan(AddPoint(TString::Format("%s=%g", v->GetName(), nextPoint)).getVal(sOpt).first)) { - if (maxTries == 0) { - ::Warning("xRooHypoSpace::findlimit", "Reached max number of point evaluations"); - } else { - ::Error("xRooHypoSpace::findlimit", "Problem evaluating %s @ %s=%g", sOpt.Data(), v->GetName(), nextPoint); - } + if (std::isnan(AddPoint(TString::Format("%s=%g", v->GetName(), nextPoint)).getVal(sOpt).first)) { + ::Error("xRooHypoSpace::findlimit", "Problem evaluating %s @ %s=%g", sOpt.Data(), v->GetName(), nextPoint); return lim; } gr.reset(); @@ -1578,7 +1604,8 @@ void xRooNLLVar::xRooHypoSpace::Draw(Option_t *opt) auto gra2 = static_cast(gra->DrawClone(sOpt.Contains("same") ? "" : "A")); gra2->SetBit(kCanDelete); if (sOpt.Contains("pcls") || sOpt.Contains("pnull")) { - gra2->GetHistogram()->SetMinimum(1e-6); + gra2->SetMinimum(1e-6); + gra2->SetMaximum(1); } if (gPad) { gPad->RedrawAxis(); @@ -1667,9 +1694,9 @@ void xRooNLLVar::xRooHypoSpace::Draw(Option_t *opt) minMax.second = std::max(minMax.second, val); } if (minMax.first < std::numeric_limits::infinity()) - out->GetHistogram()->SetMinimum(minMax.first); + out->SetMinimum(minMax.first); if (minMax.second > -std::numeric_limits::infinity()) - out->GetHistogram()->SetMaximum(minMax.second); + out->SetMaximum(minMax.second); TGraph *badPoints = nullptr; diff --git a/roofit/xroofit/src/xRooNLLVar.cxx b/roofit/xroofit/src/xRooNLLVar.cxx index c96601236e5df..d1bf059351b38 100644 --- a/roofit/xroofit/src/xRooNLLVar.cxx +++ b/roofit/xroofit/src/xRooNLLVar.cxx @@ -687,7 +687,8 @@ RooArgList xRooNLLVar::xRooFitResult::ranknp(const char *poi, bool up, bool pref auto vv = static_cast(out.at(out.size() - 1)); vv->setVal(v); vv->removeError(); - vv->removeMin();vv->removeMax();//vv->removeRange(); + vv->removeMin(); + vv->removeMax(); // vv->removeRange(); } return out; } @@ -712,7 +713,8 @@ xRooNLLVar::xRooFitResult xRooNLLVar::minimize(const std::shared_ptr(out->constPars()).addClone(RooRealVar(".pgof", "GoF p-value", pgof())); // and just main term - const_cast(out->constPars()).addClone(RooRealVar(".mainterm_pgof", "MainTerm GoF p-value", mainTermPgof())); + const_cast(out->constPars()) + .addClone(RooRealVar(".mainterm_pgof", "MainTerm GoF p-value", mainTermPgof())); } return xRooFitResult(std::make_shared(out, fPdf), std::make_shared(*this)); @@ -743,7 +745,7 @@ double xRooNLLVar::getEntryVal(size_t entry) const *std::unique_ptr(_pdf->getObservables(_data)) = *_data->get(entry); // if (auto s = dynamic_cast(_pdf.get());s) return // -_data->weight()*s->getPdf(s->indexCat().getLabel())->getLogVal(_data->get()); - return -_data->weight() * _pdf->getLogVal(_data->get()); + return (_data->weight() == 0) ? 0 : (-_data->weight() * _pdf->getLogVal(_data->get())); } std::set xRooNLLVar::binnedChannels() const @@ -1469,7 +1471,7 @@ xRooNLLVar::xValueWithError xRooNLLVar::xRooHypoPoint::getVal(const char *what) TString toyNum = sWhat(sWhat.Index("toys=") + 5, sWhat.Length()); size_t nToys = toyNum.Atoi(); size_t nToysAlt = (toyNum.Atof() - nToys) * nToys; - if (nToysAlt == 0 && !toyNum.Contains('.')) + if (nToysAlt == 0 && !toyNum.Contains('.') && !doNull) nToysAlt = nToys; if (nullToys.size() < nToys) { addNullToys(nToys - nullToys.size()); @@ -1593,7 +1595,15 @@ void xRooNLLVar::xRooHypoPoint::Print(Option_t *) const if (!std::isnan(v->getVal())) any_alt = true; } - std::cout << " , pllType: " << fPllType << std::endl; + std::cout << " , pllType: "; + switch (fPllType) { + case 0: std::cout << "tmu"; break; + case 1: std::cout << "qmu or qmutilde"; break; // should check for 'physical' to decide if is latter + case 2: std::cout << "q0"; break; + case 4: std::cout << "u0"; break; + default: std::cout << "unknown"; break; + } + std::cout << std::endl; if (fPllType == xRooFit::Asymptotics::Unknown) { std::cout << " obs ts: " << obs_ts << " +/- " << obs_ts_err << std::endl; @@ -1765,6 +1775,14 @@ std::shared_ptr xRooNLLVar::xRooHypoPoint::asimov(boo // dynamic_cast(p)->removeRange("physical"); -- can't use this as will modify shared property if (auto v = dynamic_cast(p)) { v->deleteSharedProperties(); // effectively removes all custom ranges + if (v->getVal() == 0) { + // for discovery tests, we generate asimov at mu!=0 and then evaluate the two sided + // at some value of mu. Normally we would use mu=0 but if we have a bin + // with only signal contribution (no bkg) will get asimov data in that bin + // and no prediction ... the cfit(mu=0) will never succeed on this + // so lets move to half the alt value instead (the value used to generate) + v->setVal(theFit->constPars().getRealValue(v->GetName()) * 0.5); + } } } @@ -1789,15 +1807,19 @@ xRooNLLVar::xValueWithError xRooNLLVar::xRooHypoPoint::pNull_asymp(double nSigma auto first_poi = dynamic_cast(poi().first()); if (!first_poi) return std::pair(std::numeric_limits::quiet_NaN(), 0); - auto _sigma_mu = sigma_mu(); + double lowBound = first_poi->getMin("physical"); + double hiBound = first_poi->getMax("physical"); + // don't need to calculate sigma_mu if physical boundaries at infinity, PValue doesn't depend on it + auto _sigma_mu = + (lowBound == -std::numeric_limits::infinity() && hiBound == std::numeric_limits::infinity()) + ? std::pair(0, 0) + : sigma_mu(); double nom = xRooFit::Asymptotics::PValue(fPllType, ts_asymp(nSigma).first, fNullVal(), fNullVal(), _sigma_mu.first, - first_poi->getMin("physical"), first_poi->getMax("physical")); - double up = - xRooFit::Asymptotics::PValue(fPllType, ts_asymp(nSigma).first + ts_asymp(nSigma).second, fNullVal(), fNullVal(), - _sigma_mu.first, first_poi->getMin("physical"), first_poi->getMax("physical")); - double down = - xRooFit::Asymptotics::PValue(fPllType, ts_asymp(nSigma).first - ts_asymp(nSigma).second, fNullVal(), fNullVal(), - _sigma_mu.first, first_poi->getMin("physical"), first_poi->getMax("physical")); + lowBound, hiBound); + double up = xRooFit::Asymptotics::PValue(fPllType, ts_asymp(nSigma).first + ts_asymp(nSigma).second, fNullVal(), + fNullVal(), _sigma_mu.first, lowBound, hiBound); + double down = xRooFit::Asymptotics::PValue(fPllType, ts_asymp(nSigma).first - ts_asymp(nSigma).second, fNullVal(), + fNullVal(), _sigma_mu.first, lowBound, hiBound); return std::pair(nom, std::max(std::abs(up - nom), std::abs(down - nom))); } @@ -2238,8 +2260,8 @@ xRooNLLVar::xValueWithError xRooNLLVar::xRooHypoPoint::sigma_mu(bool readOnly) } auto out = asi->pll(readOnly); - return std::pair(std::abs(fNullVal() - fAltVal()) / sqrt(out.first), - out.second * 0.5 * std::abs(fNullVal() - fAltVal()) / + return std::pair(std::abs(asi->fNullVal() - fAltVal()) / sqrt(out.first), + out.second * 0.5 * std::abs(asi->fNullVal() - fAltVal()) / (out.first * sqrt(out.first))); } @@ -2605,6 +2627,10 @@ xRooNLLVar::hypoPoint(const char *poiValues, double alt_value, const xRooFit::As AutoRestorer snap(*fFuncVars); out.nllVar = std::make_shared(*this); + // clear the underlying RooAbsReal, so that we don't accidentally alter it (e.g. setAttribute readOnly) + // and therefore alter the xRooNLLVar object we are creating this hypoPoint from + // basically ensure the hypoPoint has an independent version of the function + out.nllVar->reset(); out.fData = getData(); TStringToken pattern(poiValues, ","); @@ -2635,14 +2661,7 @@ xRooNLLVar::hypoPoint(const char *poiValues, double alt_value, const xRooFit::As if (poiNames == "") { throw std::runtime_error("No poi"); } - if (!std::isnan(alt_value)) { - std::unique_ptr thePoi(fFuncVars->selectByName(poiNames)); - for (auto b : *thePoi) { - if (!static_cast(b)->hasRange("physical")) { - static_cast(b)->setRange("physical", 0, std::numeric_limits::infinity()); - } - } - } + auto _snap = std::unique_ptr(fFuncVars->selectByAttrib("Constant", true))->snapshot(); _snap->setAttribAll("poi", false); std::unique_ptr _poi(_snap->selectByName(poiNames)); @@ -2667,11 +2686,34 @@ xRooNLLVar::hypoPoint(const char *poiValues, double alt_value, const xRooFit::As _type = xRooFit::Asymptotics::OneSidedPositive; } else { _type = xRooFit::Asymptotics::Uncapped; + // for uncapped, should check min is not at physical boundary + for (auto b : out.poi()) { + if (auto r = dynamic_cast(b)) { + if (r->hasRange("physical") && r->getMin() >= r->getMin("physical")) { + ::Info( + "xRooNLLVar::hypoPoint", + "fitting min of %s is >= physical limit (%g), but using uncapped test-statistic, so will set to " + "-max = %g", + r->GetName(), r->getMin("physical"), -r->getMax()); + r->setMin(-r->getMax()); + } + } + } } } out.fPllType = _type; + // if doing onesidedpositive with an alt value, will assume we need a physical boundary + if (!std::isnan(alt_value) && out.fPllType == xRooFit::Asymptotics::OneSidedPositive) { + std::unique_ptr thePoi(fFuncVars->selectByName(poiNames)); + for (auto b : *thePoi) { + if (!static_cast(b)->hasRange("physical")) { + static_cast(b)->setRange("physical", 0, std::numeric_limits::infinity()); + } + } + } + return out; } @@ -3055,7 +3097,8 @@ xRooNLLVar::xRooHypoSpace xRooNLLVar::hypoSpace(const char *parName, int nPoints dynamic_cast(p)->setRange("physical", 0, std::numeric_limits::infinity()); Info("xRooNLLVar::hypoSpace", "Setting physical range of %s to [0,inf]", p->GetName()); } else if (auto v = dynamic_cast(p); v->hasRange("physical")) { - v->removeMin("physical"); v->removeMax("physical");//v->removeRange("physical"); + v->removeMin("physical"); + v->removeMax("physical"); // v->removeRange("physical"); Info("xRooNLLVar::hypoSpace", "Removing physical range of %s", p->GetName()); } } @@ -3072,11 +3115,13 @@ xRooNLLVar::xRooHypoSpace xRooNLLVar::hypoSpace(const char *parName, int nPoints if (nPoints > 0) { out.AddPoints(parName, nPoints, low, high); } else { - if (!std::isnan(low) && !std::isnan(high) && !(std::isinf(low) && std::isinf(high))) { - for (auto p : out.poi()) { - dynamic_cast(p)->setRange("scan", low, high); + // if (!std::isnan(low) && !std::isnan(high) && !(std::isinf(low) && std::isinf(high))) { + for (auto p : out.poi()) { + if (auto r = dynamic_cast(p)) { + r->setRange("scan", std::isnan(low) ? r->getMin() : low, std::isnan(high) ? r->getMax() : high); } } + //} } return out; } @@ -3085,11 +3130,13 @@ xRooNLLVar::xRooHypoSpace xRooNLLVar::hypoSpace(const char *parName, int nPoints if (nPoints > 0) hs.AddPoints(parName, nPoints, low, high); else { - if (!std::isnan(low) && !std::isnan(high) && !(std::isinf(low) && std::isinf(high))) { - for (auto p : hs.poi()) { - dynamic_cast(p)->setRange("scan", low, high); + // if (!std::isnan(low) && !std::isnan(high) && !(std::isinf(low) && std::isinf(high))) { + for (auto p : hs.poi()) { + if (auto r = dynamic_cast(p)) { + r->setRange("scan", std::isnan(low) ? r->getMin() : low, std::isnan(high) ? r->getMax() : high); } } + //} } return hs; } @@ -3120,10 +3167,18 @@ xRooNLLVar::hypoSpace(const char *parName, const xRooFit::Asymptotics::PLLType & throw std::runtime_error("You must specify at least one POI for the hypoSpace"); }*/ s.fNlls[s.fPdfs.begin()->second] = std::make_shared(*this); + // clear the underlying RooAbsReal, so that we don't accidentally alter it (e.g. setAttribute readOnly) + // and therefore alter the xRooNLLVar object we are creating this hypoPoint from + // basically ensure the hypoPoint has an independent version of the function + s.fNlls[s.fPdfs.begin()->second]->reset(); s.fTestStatType = pllType; for (auto poi : s.poi()) { poi->setStringAttribute("altVal", std::isnan(alt_value) ? nullptr : TString::Format("%f", alt_value)); + // default scan range to range of poi + if (auto r = dynamic_cast(poi)) { + r->setRange("scan", r->getMin(), r->getMax()); + } } return s; @@ -3299,13 +3354,15 @@ RooStats::HypoTestResult xRooNLLVar::xRooHypoPoint::result() return out; } -std::string cling::printValue(const xRooNLLVar::xValueWithError *v) +END_XROOFIT_NAMESPACE + +std::string cling::printValue(const XROOFIT_NAMESPACE_NAME::xRooNLLVar::xValueWithError *v) { if (!v) return "xValueWithError: nullptr\n"; - return Form("%f +/- %f", v->first, v->second); + return v->__repr__(); } -std::string cling::printValue(const std::map *m) +std::string cling::printValue(const std::map *m) { if (!m) return "nullptr\n"; @@ -3315,6 +3372,4 @@ std::string cling::printValue(const std::map &comp, const "this msg)", int(noErrorPars.size()), (*noErrorPars.begin())->GetName(), (noErrorPars.size() > 1) ? ",..." : ""); if (gEnv->GetValue("XRooFit.SkipInitParErrorInference", false)) { - Warning("xRooNode", "Skipping because XRooFit.SkipInitParErrorInference=true. This is expert-only, you should fix your workspaces!"); + Warning("xRooNode", "Skipping because XRooFit.SkipInitParErrorInference=true. This is expert-only, you " + "should fix your workspaces!"); } else { // get the first top-level pdf browse(); @@ -823,16 +824,22 @@ void xRooNode::Browse(TBrowser *b) formu.ReplaceAll(TString::Format("x[%zu]", i), gv->dependents()[i].GetName()); } _name += formu; - } else if(auto pi = v->get()) { + } else if (auto pi = v->get()) { // check if all interpCodes are the same. std::set interpCodes; - for(auto& c : pi->interpolationCodes()) interpCodes.insert(c); - if(interpCodes.size()==1) { _name += TString::Format(" [InterpCode=%d]",*interpCodes.begin()); } - } else if(auto fiv = v->get()) { + for (auto &c : pi->interpolationCodes()) + interpCodes.insert(c); + if (interpCodes.size() == 1) { + _name += TString::Format(" [InterpCode=%d]", *interpCodes.begin()); + } + } else if (auto fiv = v->get()) { // check if all interpCodes are the same. std::set interpCodes; - for(auto& c : fiv->interpolationCodes()) interpCodes.insert(c==4 ? 5 : c); // in definition of FlexibleInterpVar 4 gets replaced with 5 - if(interpCodes.size()==1) { _name += TString::Format(" [InterpCode=%d]",*interpCodes.begin()); } + for (auto &c : fiv->interpolationCodes()) + interpCodes.insert(c == 4 ? 5 : c); // in definition of FlexibleInterpVar 4 gets replaced with 5 + if (interpCodes.size() == 1) { + _name += TString::Format(" [InterpCode=%d]", *interpCodes.begin()); + } } // tool tip defaults to displaying name and title, so temporarily set name to obj name if has one // and set title to the object type @@ -1357,9 +1364,10 @@ const char *xRooNode::GetIconName() const const char *xRooNode::GetNodeType() const { - if(auto rrs = get(); rrs) { + if (auto rrs = get(); rrs) { // if is BinnedLikelihood show that option - if(rrs->getAttribute("BinnedLikelihood")) return "BinnedLikelihood"; + if (rrs->getAttribute("BinnedLikelihood")) + return "BinnedLikelihood"; } if (auto o = get(); o && fParent && (fParent->get() || fParent->get())) { if (o->InheritsFrom("RooStats::HistFactory::FlexibleInterpVar")) @@ -1536,7 +1544,7 @@ xRooNode xRooNode::Remove(const xRooNode &child) arg = p2->components().find(child.GetName()); if (!arg) throw std::runtime_error(TString::Format("Cannot find %s in %s", child.GetName(), fParent->GetName())); - // remove server ... doesn't seem to trigger removal from proxy + // remove server ... doesn't seem to trigger removal from proxy #if ROOT_VERSION_CODE < ROOT_VERSION(6, 27, 00) p2->_compRSet.remove(*arg); #else @@ -1916,11 +1924,11 @@ xRooNode xRooNode::Add(const xRooNode &child, Option_t *opt) return *this; } auto _arg = child.get(); - if(auto _ds = dynamic_cast(p); _arg && _ds) { + if (auto _ds = dynamic_cast(p); _arg && _ds) { // can add var or function of existing obs to dataset as a column _ds->addColumn(*_arg); _arg->setAttribute("obs"); - return xRooNode(*_arg,*this); + return xRooNode(*_arg, *this); } auto _h = child.get(); if (!_h) { @@ -1992,7 +2000,7 @@ xRooNode xRooNode::Add(const xRooNode &child, Option_t *opt) if (auto p = get(); p) { auto cc = child.fComp; - //bool isConverted = (cc != child.convertForAcquisition(*this, sOpt)); + // bool isConverted = (cc != child.convertForAcquisition(*this, sOpt)); child.convertForAcquisition(*this, sOpt); if ((child.get() || (!child.fComp && getObject(child.GetName())))) { auto out = (child.fComp) ? acquire(child.fComp) : getObject(child.GetName()); @@ -2810,16 +2818,22 @@ void xRooNode::Print(Option_t *opt) const formu.ReplaceAll(TString::Format("x[%zu]", i), gv->dependents()[i].GetName()); } _suffix += formu; - } else if(auto pi = get()) { + } else if (auto pi = get()) { // check if all interpCodes are the same. Will include in the NodeType std::set interpCodes; - for(auto& c : pi->interpolationCodes()) interpCodes.insert(c); - if(interpCodes.size()==1) { _suffix += TString::Format(" [InterpCode=%d]",*interpCodes.begin()); } - } else if(auto fiv = get()) { + for (auto &c : pi->interpolationCodes()) + interpCodes.insert(c); + if (interpCodes.size() == 1) { + _suffix += TString::Format(" [InterpCode=%d]", *interpCodes.begin()); + } + } else if (auto fiv = get()) { // check if all interpCodes are the same. std::set interpCodes; - for(auto& c : fiv->interpolationCodes()) interpCodes.insert(c==4 ? 5 : c); // in definition of FlexibleInterpVar 4 gets replaced with 5 - if(interpCodes.size()==1) { _suffix += TString::Format(" [InterpCode=%d]",*interpCodes.begin()); } + for (auto &c : fiv->interpolationCodes()) + interpCodes.insert(c == 4 ? 5 : c); // in definition of FlexibleInterpVar 4 gets replaced with 5 + if (interpCodes.size() == 1) { + _suffix += TString::Format(" [InterpCode=%d]", *interpCodes.begin()); + } } std::cout << get()->ClassName() << "::" << get()->GetName() << _suffix.Data() << std::endl; } @@ -2880,16 +2894,22 @@ void xRooNode::Print(Option_t *opt) const formu.ReplaceAll(TString::Format("x[%zu]", j), gv->dependents()[j].GetName()); } _suffix += formu; - } else if(auto pi = k->get()) { + } else if (auto pi = k->get()) { // check if all interpCodes are the same. Will include in the NodeType std::set interpCodes; - for(auto& c : pi->interpolationCodes()) interpCodes.insert(c); - if(interpCodes.size()==1) { _suffix += TString::Format(" [InterpCode=%d]",*interpCodes.begin()); } - } else if(auto fiv = k->get()) { + for (auto &c : pi->interpolationCodes()) + interpCodes.insert(c); + if (interpCodes.size() == 1) { + _suffix += TString::Format(" [InterpCode=%d]", *interpCodes.begin()); + } + } else if (auto fiv = k->get()) { // check if all interpCodes are the same. std::set interpCodes; - for(auto& c : fiv->interpolationCodes()) interpCodes.insert(c==4 ? 5 : c); // in definition of FlexibleInterpVar 4 gets replaced with 5 - if(interpCodes.size()==1) { _suffix += TString::Format(" [InterpCode=%d]",*interpCodes.begin()); } + for (auto &c : fiv->interpolationCodes()) + interpCodes.insert(c == 4 ? 5 : c); // in definition of FlexibleInterpVar 4 gets replaced with 5 + if (interpCodes.size() == 1) { + _suffix += TString::Format(" [InterpCode=%d]", *interpCodes.begin()); + } } std::cout << k->get()->ClassName() << "::" << k->get()->GetName() << _suffix.Data() << std::endl; } @@ -4234,7 +4254,8 @@ void xRooNode::_fit_(const char *constParValues, const char *options) : gClient->GetRoot(); TString gofResult = ""; if (_nll.fOpts->find("GoF")) { - gofResult = TString::Format("GoF p-value = %g (mainTerm = %g)\n", fr->constPars().getRealValue(".pgof"),fr->constPars().getRealValue(".mainterm_pgof")); + gofResult = TString::Format("GoF p-value = %g (mainTerm = %g)\n", fr->constPars().getRealValue(".pgof"), + fr->constPars().getRealValue(".mainterm_pgof")); } if (fr->status() != 0) { new TGMsgBox(gClient->GetRoot(), w, "Fit Finished with Bad Status Code", @@ -4943,7 +4964,7 @@ bool xRooNode::SetBinError(int bin, double value) TString origName = (f->getStringAttribute("origName")) ? f->getStringAttribute("origName") : GetName(); rrv->setStringAttribute(Form("sumw2_%s", origName.Data()), TString::Format("%f", pow(value, 2))); auto bin_pars = f->dataHist().get(bin - 1); - auto _binContent = f->dataHist().weight(bin-1); + auto _binContent = f->dataHist().weight(bin - 1); if (f->getAttribute("density")) { _binContent *= f->dataHist().binVolume(*bin_pars); } @@ -5299,9 +5320,10 @@ std::shared_ptr xRooNode::convertForAcquisition(xRooNode &acquirer, con fComp = _f; return _f; - } else if (!get() && (sName.BeginsWith("factory:")||sName.Contains("::")) && acquirer.ws()) { + } else if (!get() && (sName.BeginsWith("factory:") || sName.Contains("::")) && acquirer.ws()) { TString s(sName); - if(sName.BeginsWith("factory:")) s = TString(s(8, s.Length())); + if (sName.BeginsWith("factory:")) + s = TString(s(8, s.Length())); fComp.reset(acquirer.ws()->factory(s), [](TObject *) {}); if (fComp) { const_cast(this)->TNamed::SetName(fComp->GetName()); @@ -6316,14 +6338,11 @@ xRooNode xRooNode::vars() const } } } else if (auto w = get(); w) { - for (auto a : w->allVars()) { - out.emplace_back(std::make_shared(*a, *this)); - out.get()->add(*a); - } - // add all cats as well - for (auto a : w->allCats()) { - out.emplace_back(std::make_shared(*a, *this)); - out.get()->add(*a); + for (auto a : w->components()) { + if (a->InheritsFrom(RooRealVar::Class()) || a->InheritsFrom(RooCategory::Class()) || + a->InheritsFrom(RooConstVar::Class())) { + out.emplace_back(std::make_shared(*a, *this)); + } } } return out; @@ -6930,6 +6949,33 @@ xRooNode xRooNode::datasets() const return out; } +xRooNode xRooNode::parents() const +{ + xRooNode out(".parents", nullptr, *this); + if (auto a = get()) { + for (auto c : a->clients()) { + out.push_back(std::make_shared(*c, *this)); + } + } + return out; +} + +xRooNode xRooNode::args() const +{ + if (auto w = get()) { + xRooNode out(".args", w->components(), *this); + out.browse(); // populate + return out; + } else if (auto a = get()) { + xRooNode out(".args", std::make_shared(), *this); + out.get()->setName((GetPath() + ".args").c_str()); + a->treeNodeServerList(out.get()); + out.browse(); // populate + return out; + } + return nullptr; +} + std::shared_ptr xRooNode::getBrowsable(const char *name) const { for (auto b : fBrowsables) { @@ -7021,8 +7067,21 @@ TGraph *xRooNode::BuildGraph(RooAbsLValue *v, bool includeZeros, TVirtualPad *fr dataGraph->SetTitle(TString::Format("%s;%s;Events", dataGraph->GetTitle(), theHist->GetXaxis()->GetTitle())); *static_cast(dataGraph) = *static_cast(theHist); *static_cast(dataGraph) = *static_cast(theHist); + + // default style based on if generated or not + + if (auto w = theData->weightVar(); w && w->getStringAttribute("fitResult")) { + // is generated + dataGraph->SetLineColor(kGreen + 2); + if (w->getAttribute("expected")) { + // is asimov + dataGraph->SetLineColor(kBlue); + } + } else { + dataGraph->SetLineColor(kBlack); + } dataGraph->SetMarkerStyle(20); - dataGraph->SetLineColor(kBlack); + dataGraph->SetMarkerColor(dataGraph->GetLineColor()); dataGraph->SetMarkerSize(gStyle->GetMarkerSize()); auto _obs = obs(); @@ -7902,7 +7961,10 @@ xRooNode xRooNode::reduced(const std::string &_range, bool invert) const _cat.setLabel(cName); bool matchAny = false; for (auto &p : patterns) { - if (cName.Contains(TRegexp(p, true))) { + TString pNoCatName(p); + if (pNoCatName.Contains('=')) + pNoCatName = pNoCatName(pNoCatName.Index('=') + 1, pNoCatName.Length()); + if (cName.Contains(TRegexp(p, true)) || cName.Contains(TRegexp(pNoCatName, true))) { matchAny = true; break; } @@ -8006,6 +8068,26 @@ xRooNode xRooNode::reduced(const std::string &_range, bool invert) const return get() ? xRooNode(std::make_shared(), fParent) : *this; } +xRooNode xRooNode::reduced(const std::function selector) const +{ + if (empty()) { + const_cast(*this).browse(); + } + // build a list of children to keep + std::string noName = "___"; // assume this is never a name + std::string childNames; + for (auto &c : *this) { + if (selector(*c)) { + if (!childNames.empty()) + childNames += ","; + childNames += c->GetName(); + } + } + if (childNames.empty()) + childNames = noName; // if childNames was blank it would return the full list; + return reduced(childNames); // calls main method above ... this will ensure we construct a reduced version of ourself +} + // xRooNode xRooNode::generate(bool expected) const { // // auto fr = fitResult(); @@ -8099,26 +8181,18 @@ class xRooProjectedPdf : public RooProjectedPdf { } }; -double new_getPropagatedError(const RooAbsReal &f, const RooFitResult &fr, const RooArgSet &nset = {}, - RooArgList **pars = nullptr, bool asymHi = false, bool asymLo = false) +std::vector new_getPropagatedErrors(const RooAbsReal &f, const RooFitResult &fr, const RooArgSet &nset, + RooArgList **pars, bool asymHi, bool asymLo, const std::vector &bins, + const std::function &setBin) { - // Calling getParameters() might be costly, but necessary to get the right - // parameters in the RooAbsReal. The RooFitResult only stores snapshots. - - // handle simple case that function is a RooRealVar - if (auto rrv = dynamic_cast(&f); rrv) { - if (auto frrrv = dynamic_cast(fr.floatParsFinal().find(*rrv)); frrrv) { - rrv = frrrv; // use value from fit result - } - if (asymHi) { - return rrv->getErrorHi(); - } else if (asymLo) { - return rrv->getErrorLo(); - } else { - return rrv->getError(); - } - } + const size_t nBins = bins.size(); + std::vector out(nBins, 0.); + if (nBins == 0) + return out; + // Build the list of parameters f depends on that have a non-zero error. + // The result is cached in *pars so it is only computed once across repeated + // calls (e.g. the separate Hi and Lo passes of an asymmetric error band). RooArgList *_pars = (pars) ? *pars : nullptr; if (!_pars) { @@ -8155,83 +8229,119 @@ double new_getPropagatedError(const RooAbsReal &f, const RooFitResult &fr, const } } - // Make std::vector of variations - TVectorD F(_pars->size()); - - // Create std::vector of plus,minus variations for each parameter - TMatrixDSym V(_pars->size() == fr.floatParsFinal().size() ? fr.covarianceMatrix() - : fr.reducedCovarianceMatrix(*_pars)); + const size_t nPars = _pars->size(); + // Covariance and correlation matrices depend only on the fit result and the + // parameter set, so build them (including the matrix reduction) once for all + // bins rather than once per bin. // TODO: if _pars includes pars not in fr, need to extend matrix with uncorrelated errors of those pars + // (as built above, _pars is a subset of fr.floatParsFinal(), so such pars are currently dropped silently). + TMatrixDSym V(nPars == fr.floatParsFinal().size() ? fr.covarianceMatrix() : fr.reducedCovarianceMatrix(*_pars)); - double nomVal = f.getVal(nset); + TMatrixDSym C(nPars); + for (size_t i = 0; i < nPars; i++) { + for (size_t j = i; j < nPars; j++) { + C(i, j) = V(i, j) / std::sqrt(V(i, i) * V(j, j)); + C(j, i) = C(i, j); + } + } + + // Nominal value of each bin (only needed to symmetrize asymmetric errors). + std::vector nomVals; + if (asymHi || asymLo) { + nomVals.resize(nBins); + for (size_t k = 0; k < nBins; k++) { + setBin(bins[k]); + nomVals[k] = f.getVal(nset); + } + } - for (std::size_t ivar = 0; ivar < _pars->size(); ivar++) { + // F[k] holds, for bin k, the vector of (signed) variations w.r.t. each + // parameter - i.e. the column of the Jacobian times the parameter error. + std::vector F(nBins, TVectorD(nPars)); + + std::vector plusVar(nBins), minusVar(nBins); + + for (size_t ivar = 0; ivar < nPars; ivar++) { auto &rrv = static_cast((*_pars)[ivar]); auto *frrrv = static_cast(fr.floatParsFinal().find(rrv)); double cenVal = rrv.getVal(); - double plusVar, minusVar, errVal; if (asymHi || asymLo) { - errVal = frrrv->getErrorHi(); - rrv.setVal(cenVal + errVal); - plusVar = f.getVal(nset); - errVal = frrrv->getErrorLo(); - rrv.setVal(cenVal + errVal); - minusVar = f.getVal(nset); - if (asymHi) { - // pick the one that moved result 'up' most - plusVar = std::max(plusVar, minusVar); - minusVar = 2 * nomVal - plusVar; // symmetrizes - } else { - // pick the one that moved result 'down' most - minusVar = std::min(plusVar, minusVar); - plusVar = 2 * nomVal - minusVar; // symmetrizes + double errValHi = frrrv->getErrorHi(); + rrv.setVal(errValHi > 0 ? std::min(cenVal + errValHi, rrv.getMax()) + : std::max(cenVal + errValHi, rrv.getMin())); + for (size_t k = 0; k < nBins; k++) { + setBin(bins[k]); + plusVar[k] = f.getVal(nset); + } + double errValLo = frrrv->getErrorLo(); + rrv.setVal(errValLo > 0 ? std::min(cenVal + errValLo, rrv.getMax()) + : std::max(cenVal + errValLo, rrv.getMin())); + for (size_t k = 0; k < nBins; k++) { + setBin(bins[k]); + minusVar[k] = f.getVal(nset); + } + rrv.setVal(cenVal); + for (size_t k = 0; k < nBins; k++) { + double hi = plusVar[k]; + double lo = minusVar[k]; + if (asymHi) { + // pick the one that moved result 'up' most + hi = std::max(plusVar[k], minusVar[k]); + lo = 2 * nomVals[k] - hi; // symmetrizes + } else { + // pick the one that moved result 'down' most + lo = std::min(plusVar[k], minusVar[k]); + hi = 2 * nomVals[k] - lo; // symmetrizes + } + F[k][ivar] = (hi - lo) * 0.5; } } else { - errVal = sqrt(V(ivar, ivar)); + double errVal = sqrt(V(ivar, ivar)); // Make Plus variation - rrv.setVal(cenVal + errVal); - plusVar = f.getVal(nset); + rrv.setVal(std::min(cenVal + errVal, rrv.getMax())); + for (size_t k = 0; k < nBins; k++) { + setBin(bins[k]); + plusVar[k] = f.getVal(nset); + } // Make Minus variation - rrv.setVal(cenVal - errVal); - minusVar = f.getVal(nset); + rrv.setVal(std::max(cenVal - errVal, rrv.getMin())); + for (size_t k = 0; k < nBins; k++) { + setBin(bins[k]); + minusVar[k] = f.getVal(nset); + } + rrv.setVal(cenVal); + for (size_t k = 0; k < nBins; k++) { + F[k][ivar] = (plusVar[k] - minusVar[k]) * 0.5; + } } - F[ivar] = (plusVar - minusVar) * 0.5; - rrv.setVal(cenVal); } - // Re-evaluate this RooAbsReal with the central parameters just to be - // extra-safe that a call to `getPropagatedError()` doesn't change any state. - // It should not be necessary because thanks to the dirty flag propagation - // the RooAbsReal is re-evaluated anyway the next time getVal() is called. - // Still there are imaginable corner cases where it would not be triggered, - // for example if the user changes the RooFit operation more after the error - // propagation. + // Re-evaluate f with the central parameters just to be extra-safe that this + // call doesn't leave any state changed. It should not be necessary because + // thanks to the dirty flag propagation f is re-evaluated anyway the next + // time getVal() is called, but there are imaginable corner cases where that + // would not be triggered (e.g. if the caller changes the RooFit operation + // mode after the error propagation). + setBin(bins[nBins - 1]); f.getVal(nset); - TMatrixDSym C(_pars->size()); - std::vector errVec(_pars->size()); - for (std::size_t i = 0; i < _pars->size(); i++) { - errVec[i] = std::sqrt(V(i, i)); - for (std::size_t j = i; j < _pars->size(); j++) { - C(i, j) = V(i, j) / std::sqrt(V(i, i) * V(j, j)); - C(j, i) = C(i, j); - } + // Calculate error in linear approximation from variations and correlation + // coefficients - this is pure arithmetic (no function evaluations). + for (size_t k = 0; k < nBins; k++) { + out[k] = std::sqrt(F[k] * (C * F[k])); } - // Calculate error in linear approximation from variations and correlation coefficient - double sum = F * (C * F); - if (!pars) { delete _pars; } else { *pars = _pars; } - return sqrt(sum); + return out; } class PdfWrapper : public RooAbsPdf { @@ -8256,7 +8366,7 @@ class PdfWrapper : public RooAbsPdf { } fExpectedEventsMode = expEvMode; } - ~PdfWrapper() override{}; + ~PdfWrapper() override {}; PdfWrapper(const PdfWrapper &other, const char *name = nullptr) : RooAbsPdf(other, name), fFunc("func", this, other.fFunc), @@ -8621,7 +8731,6 @@ TH1 *xRooNode::BuildHistogram(RooAbsLValue *v, bool empty, bool errors, int binS TObject *vv = rar; - TH1 *h = nullptr; if (!v) { if (binStart != -1 || binEnd != -1) { // allow v to stay nullptr if doing integral (binStart=binEnd=-1) @@ -8650,8 +8759,6 @@ TH1 *xRooNode::BuildHistogram(RooAbsLValue *v, bool empty, bool errors, int binS TDirectory::TContext ctx{nullptr}; // No self-registration to directories h = static_cast(templateHist->Clone(rar->GetName())); } - if (h->GetListOfFunctions()) - h->GetListOfFunctions()->Clear(); h->SetTitle(rar->GetTitle()); h->Reset(); } else if (x) { @@ -8872,13 +8979,17 @@ TH1 *xRooNode::BuildHistogram(RooAbsLValue *v, bool empty, bool errors, int binS for (auto o : _obs) { if (auto rr = o->get(); rr && rr->hasRange("coordRange")) { - rr->removeMin();rr->removeMax();//rr->removeRange("coordRange"); // doesn't actually remove, just sets to -inf->+inf + rr->removeMin(); + rr->removeMax(); // rr->removeRange("coordRange"); // doesn't actually remove, just sets to + // -inf->+inf rr->setStringAttribute("coordRange", nullptr); // removes the attribute } } // probably should also remove any range on the x-axis variable too, if there is one if (auto rr = dynamic_cast(v); rr && rr->hasRange("coordRange")) { - rr->removeMin();rr->removeMax();//rr->removeRange("coordRange"); // doesn't actually remove, just sets to -inf->+inf + rr->removeMin(); + rr->removeMax(); // rr->removeRange("coordRange"); // doesn't actually remove, just sets to + // -inf->+inf rr->setStringAttribute("coordRange", nullptr); // removes the attribute } coords(); // loads current coordinates and populates coordRange, if any @@ -9044,6 +9155,7 @@ TH1 *xRooNode::BuildHistogram(RooAbsLValue *v, bool empty, bool errors, int binS empty = false; // must not be empty b.c. calculation of error relies on knowing nominal (see after loop) } + std::vector errorBins; // bins (in-range) for which to compute the propagated error band, filled below for (int toy = 0; toy < (nErrorToys + 1); toy++) { TH1 *main_h = h; @@ -9101,64 +9213,9 @@ TH1 *xRooNode::BuildHistogram(RooAbsLValue *v, bool empty, bool errors, int binS if (errors) { static_cast(h->FindObject("nominal"))->SetBinContent(i, r); // transfer nominal to nominal hist - double res; - bool doAsym = (errorsHi && errorsLo); - if (doAsym) { - errorsHi = false; - } - if (p) { - // std::cout << "computing error of :" << h->GetBinCenter(i) << std::endl; - // //fr->floatParsFinal().Print(); fr->covarianceMatrix().Print(); - // res = PdfWrapper((oldrar) ? *rar : *p, _coefs.get(), !v, oldrar ? p : nullptr) - // .getSimplePropagatedError(*fr, normSet); - res = new_getPropagatedError( - PdfWrapper((oldrar) ? *rar : *p, _coefs.get(), !v, oldrar ? p : nullptr), *fr, normSet, - &errorPars, errorsHi, errorsLo); -#if ROOT_VERSION_CODE < ROOT_VERSION(6, 27, 00) - // improved normSet invalidity checking, so assuming no longer need this in 6.28 onwards - p->_normSet = nullptr; -#endif - } else { - // res = RooProduct("errorEval", "errorEval", - // RooArgList(*rar, !_coefs.get() ? RooFit::RooConst(1) : - // *_coefs.get())) - // .getPropagatedError( - // *fr /*, normSet*/); // should be no need to pass a normSet to a non-pdf (but - // not verified this) - res = new_getPropagatedError( - RooProduct("errorEval", "errorEval", - RooArgList(*rar, !_coefs.get() ? RooFit::RooConst(1) : *_coefs.get())), - *fr, {}, &errorPars, errorsHi, - errorsLo); // should be no need to pass a normSet to a non-pdf (but not verified this) - // especially important not to pass in the case we are evaluated RooRealSumPdf as a function! otherwise - // error will be wrong - } - if (needBinWidth) { - res *= h->GetBinWidth(i); - } - h->SetBinError(i, res); - if (doAsym) { - // compute Hi error - errorsHi = true; - errorsLo = false; - if (p) { - res = new_getPropagatedError( - PdfWrapper((oldrar) ? *rar : *p, _coefs.get(), !v, oldrar ? p : nullptr), *fr, normSet, - &errorPars, errorsHi, errorsLo); - } else { - res = new_getPropagatedError( - RooProduct("errorEval", "errorEval", - RooArgList(*rar, !_coefs.get() ? RooFit::RooConst(1) : *_coefs.get())), - *fr, {}, &errorPars, errorsHi, errorsLo); - } - if (needBinWidth) { - res *= h->GetBinWidth(i); - } - errorsLo = true; - // lowVal = content - error, highVal = content + res - // => band/2 = (res+error)/2 and band-mid = (2*content+res-error)/2 - h->SetBinContent(i, h->GetBinContent(i) + (res - h->GetBinError(i)) * 0.5); - h->SetBinError(i, (res + h->GetBinError(i)) * 0.5); + // Record for which bins to compute the error later. + if (toy == 0) { + errorBins.push_back(i); } } timeIt.Stop(); @@ -9190,6 +9247,69 @@ TH1 *xRooNode::BuildHistogram(RooAbsLValue *v, bool empty, bool errors, int binS h = main_h; } } + + if (errors && !errorBins.empty()) { + // Batched linear error propagation over all (in-range) bins recorded above. By varying each parameter only + // once and sweeping the observable inside that variation, the normalization / expected-events integrals are + // reused across bins instead of being recomputed for every bin (see new_getPropagatedErrors). + std::unique_ptr errFunc; + RooArgSet emptyNormSet; + if (p) { + errFunc = + std::make_unique((oldrar) ? *rar : *p, _coefs.get(), !v, oldrar ? p : nullptr); + } else { + // especially important not to pass a normSet in the case we are evaluating a RooRealSumPdf as a function! + // otherwise the error will be wrong + errFunc = std::make_unique( + "errorEval", "errorEval", + RooArgList(*rar, !_coefs.get() ? RooFit::RooConst(1) : *_coefs.get())); + } + const RooArgSet &errNormSet = p ? normSet : emptyNormSet; + + auto setBin = [&](int i) { + if (x) { + x->setVal(h->GetBinCenter(i)); + } else if (cat) { + cat->setLabel(h->GetXaxis()->GetBinLabel(i)); + } else if (v) { + v->setBin(i - 1); + } + }; + + bool doAsym = (errorsHi && errorsLo); + // For asymmetric bands we need both the Lo and Hi propagated errors; otherwise a single (possibly one-sided) + // pass. errorPars is reused across both passes so the parameter list is only built once. + std::vector errLo = new_getPropagatedErrors( + *errFunc, *fr, errNormSet, &errorPars, doAsym ? false : errorsHi, doAsym ? true : errorsLo, errorBins, setBin); + std::vector errHi; + if (doAsym) { + errHi = new_getPropagatedErrors(*errFunc, *fr, errNormSet, &errorPars, true, false, errorBins, setBin); + } +#if ROOT_VERSION_CODE < ROOT_VERSION(6, 27, 00) + if (p) { + // improved normSet invalidity checking, so assuming no longer need this in 6.28 onwards + p->_normSet = nullptr; + } +#endif + + for (size_t k = 0; k < errorBins.size(); k++) { + int i = errorBins[k]; + double resLo = errLo[k]; + if (needBinWidth) + resLo *= h->GetBinWidth(i); + h->SetBinError(i, resLo); + if (doAsym) { + double resHi = errHi[k]; + if (needBinWidth) + resHi *= h->GetBinWidth(i); + // lowVal = content - resLo, highVal = content + resHi + // => band/2 = (resHi+resLo)/2 and band-mid = content + (resHi-resLo)/2 + h->SetBinContent(i, h->GetBinContent(i) + (resHi - resLo) * 0.5); + h->SetBinError(i, (resHi + resLo) * 0.5); + } + } + } + if (gOldHandlerr) { signal(SIGINT, gOldHandlerr); gOldHandlerr = nullptr; @@ -9294,7 +9414,7 @@ TH1 *xRooNode::BuildHistogram(RooAbsLValue *v, bool empty, bool errors, int binS bool titleMatchName = true; std::map histGroups; std::vector hhs; - std::set> ordered_hhs; + std::set> ordered_hhs; std::set histsWithBadTitles; // these histograms will have their titles autoFormatted // support for CMS model case where has single component containing many coeffs @@ -9323,7 +9443,7 @@ TH1 *xRooNode::BuildHistogram(RooAbsLValue *v, bool empty, bool errors, int binS // seems I have to remake the function each time, as haven't figured out what cache needs clearing? zero.setAttribute(Form("ORIGNAME:%s", c->GetName())); // used in redirectServers to say what this replaces - forig->redirectServers(RooArgSet(zero), false, true); // each time will replace one additional coef + forig->redirectServers(RooArgSet(zero), false, true); // each time will replace one additional coef std::unique_ptr f(dynamic_cast(forig->Clone("tmpCopy"))); // zero.setAttribute(Form("ORIGNAME:%s",c->GetName()),false); (commented out so that on next iteration @@ -9336,8 +9456,8 @@ TH1 *xRooNode::BuildHistogram(RooAbsLValue *v, bool empty, bool errors, int binS hh->SetTitle(c->GetName()); // ensure all hists has titles // special case for CMS ... if find "_proc_" in the title, take whatever is after that auto idx = TString(hh->GetTitle()).Index("_proc_"); - if(idx>=0) { - hh->SetTitle(TString(TString(hh->GetTitle())(idx+6,strlen(hh->GetTitle())))); + if (idx >= 0) { + hh->SetTitle(TString(TString(hh->GetTitle())(idx + 6, strlen(hh->GetTitle())))); } histsWithBadTitles.insert(hh); } else if (strcmp(hh->GetName(), hh->GetTitle()) == 0) { @@ -9350,7 +9470,7 @@ TH1 *xRooNode::BuildHistogram(RooAbsLValue *v, bool empty, bool errors, int binS hh->Scale(-1.); // remove the errors ... the above lines will have introduced errors hh->TH1::Reset("ICE"); // calling the base class method explicitly will only clear errors - ordered_hhs.insert(std::pair(ordered_hhs.size(),hh)); + ordered_hhs.insert(std::pair(ordered_hhs.size(), hh)); prevHist = nextHist; } } else if (get()) { @@ -9379,7 +9499,7 @@ TH1 *xRooNode::BuildHistogram(RooAbsLValue *v, bool empty, bool errors, int binS titleMatchName &= (TString(samp->GetName()) == hh->GetTitle() || TString(hh->GetTitle()).BeginsWith(TString(samp->GetName()) + "_")); hh->SetBinContent(hh->GetXaxis()->FindFixBin(chanName), samp->GetContent()); - ordered_hhs.insert(std::pair(ordered_hhs.size(),hh)); + ordered_hhs.insert(std::pair(ordered_hhs.size(), hh)); } } } else { @@ -9390,7 +9510,11 @@ TH1 *xRooNode::BuildHistogram(RooAbsLValue *v, bool empty, bool errors, int binS hh->SetName(samp->GetName()); if (sf) hh->Scale(sf->getVal()); - ordered_hhs.insert(std::pair((samp->get() && samp->get()->getStringAttribute("StackOrder")) ? TString(samp->get()->getStringAttribute("StackOrder")).Atoi() : ordered_hhs.size(),hh)); + ordered_hhs.insert( + std::pair((samp->get() && samp->get()->getStringAttribute("StackOrder")) + ? TString(samp->get()->getStringAttribute("StackOrder")).Atoi() + : ordered_hhs.size(), + hh)); if (strlen(hh->GetTitle()) == 0) { hh->SetTitle(samp->GetName()); // ensure all hists has titles histsWithBadTitles.insert(hh); @@ -9403,7 +9527,7 @@ TH1 *xRooNode::BuildHistogram(RooAbsLValue *v, bool empty, bool errors, int binS } // pull histograms in their order - for(auto& [_,hh] : ordered_hhs) { + for (auto &[_, hh] : ordered_hhs) { hhs.push_back(hh); } @@ -10917,7 +11041,6 @@ void xRooNode::Draw(Option_t *opt) histCopy->SetBit(kCanDelete); auto _axis = (doHorizontal ? histCopy->GetYaxis() : histCopy->GetXaxis()); - graph->GetHistogram()->GetXaxis()->Set(std::max(graph->GetN(), 1), -0.5, std::max(graph->GetN(), 1) - 0.5); for (int ii = 1; ii <= _axis->GetNbins(); ii++) { graph->GetHistogram()->GetXaxis()->SetBinLabel(ii, _axis->GetBinLabel(ii)); @@ -11238,22 +11361,32 @@ void xRooNode::Draw(Option_t *opt) // drawing dataset associated to a simultaneous means must find subpads with variation names // may not have subpads if drawning a "Yield" plot ... bool doneDraw = false; - for (auto c : s->bins()) { - auto _pad = dynamic_cast(gPad->GetPrimitive(c->GetName())); + // in the case of hybrid datasets, the parentPdf will be a reducedPdf ... + // but if the dataset has been added to, then there may be additional entries + // in other channels ... so loop over all labels of the categorical + // and for ones we don't have a channel for, just draw directly on + + for (auto [catName, catVal] : s->get()->indexCat()) { + auto _pad = dynamic_cast(gPad->GetPrimitive( + TString::Format("%s=%s", s->get()->indexCat().GetName(), catName.c_str()))); if (!_pad) continue; // channel was hidden? // attach as a child before calling datasets(), so that if this dataset is external to workspace it is // included still attaching the dataset ensures dataset reduction for the channel is applied - c->push_back(std::make_shared(*this)); - auto ds = c->datasets().find(GetName()); - c->resize(c->size() - 1); // remove the child we attached - if (!ds) { - std::cout << " no ds " << GetName() << " - this should never happen!" << std::endl; - continue; - } auto tmp = gPad; _pad->cd(); - ds->Draw(opt); + if (auto c = s->bins().find(catName)) { + c->push_back(std::make_shared(*this)); + auto ds = c->datasets().find(GetName()); + c->resize(c->size() - 1); // remove the child we attached + if (!ds) { + std::cout << " no ds " << GetName() << " - this should never happen!" << std::endl; + continue; + } + ds->Draw(opt); + } else { + Draw(opt); + } doneDraw = true; tmp->cd(); } @@ -11649,7 +11782,7 @@ void xRooNode::Draw(Option_t *opt) if (!hasSame) clearPad(); - if(rar->getAttribute("Logy")) { + if (rar->getAttribute("Logy")) { gPad->SetLogy(1); } @@ -12655,7 +12788,9 @@ xRooNode::GetBinErrors(int binStart, int binEnd, const xRooNode &_fr, int nToys, // return out; } -std::string cling::printValue(const xRooNode *v) +END_XROOFIT_NAMESPACE + +std::string cling::printValue(const XROOFIT_NAMESPACE_NAME::xRooNode *v) { if (!v) return "nullptr\n"; @@ -12687,5 +12822,3 @@ std::string cling::printValue(const xRooNode *v) return out; } - -END_XROOFIT_NAMESPACE From d493eb169a82eb3b6d3991c8e2a10f33ee9a5982 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Fri, 19 Jun 2026 09:53:08 +0200 Subject: [PATCH 2/5] [RF] Fix various memory leaks caused by not deleting return values Fix various memory leaks from cases where functions that return an owned return value were used without deleting the return value. This is addressed by using `std::unique_ptr`. (cherry picked from commit 2b876bb906bbf2099ff2397b2edaeaabd8b60046) --- roofit/histfactory/src/HistFactoryModelUtils.cxx | 2 +- roofit/roofit/test/testFitPerf.cxx | 4 ++-- roofit/roofit/test/vectorisedPDFs/VectorisedPDFTests.cxx | 2 +- roofit/roofitcore/src/RooAbsReal.cxx | 6 +++--- roofit/roofitcore/src/RooRealMPFE.cxx | 3 ++- roofit/roofitcore/src/RooVectorDataStore.cxx | 5 +++-- roofit/roofitcore/src/TestStatistics/RooAbsL.cxx | 2 +- roofit/roostats/src/AsymptoticCalculator.cxx | 4 ++-- roofit/roostats/src/HypoTestInverter.cxx | 3 ++- roofit/roostats/src/LikelihoodIntervalPlot.cxx | 6 ++++-- roofit/roostats/src/PdfProposal.cxx | 5 +++-- roofit/roostats/src/SPlot.cxx | 3 ++- roofit/roostats/src/SimpleLikelihoodRatioTestStat.cxx | 6 ++++-- roofit/roostats/src/ToyMCImportanceSampler.cxx | 5 +---- roofit/xroofit/src/xRooNode.cxx | 6 +++--- 15 files changed, 34 insertions(+), 28 deletions(-) diff --git a/roofit/histfactory/src/HistFactoryModelUtils.cxx b/roofit/histfactory/src/HistFactoryModelUtils.cxx index e0a163525d2e8..5d20d1b80db67 100644 --- a/roofit/histfactory/src/HistFactoryModelUtils.cxx +++ b/roofit/histfactory/src/HistFactoryModelUtils.cxx @@ -56,7 +56,7 @@ namespace HistFactory{ if( ! FoundSumPdf ) { if(verbose) { std::cout << "Failed to find RooRealSumPdf for channel: " << sim_channel->GetName() << std::endl; - sim_channel->getComponents()->Print("V"); + components->Print("V"); } sum_pdf=nullptr; //throw std::runtime_error("Failed to find RooRealSumPdf for channel"); diff --git a/roofit/roofit/test/testFitPerf.cxx b/roofit/roofit/test/testFitPerf.cxx index da43a1df5d587..0228bfa8a98ec 100644 --- a/roofit/roofit/test/testFitPerf.cxx +++ b/roofit/roofit/test/testFitPerf.cxx @@ -551,7 +551,7 @@ int FitUsingRooFit(TTree *tree, TF1 *func) int level = 3; std::cout << "num entries = " << data.numEntries() << std::endl; bool save = true; - (pdf.getVariables())->Print("v"); // print the parameters + std::unique_ptr{pdf.getVariables()}->Print("v"); // print the parameters #else int level = -1; bool save = false; @@ -637,7 +637,7 @@ int FitUsingRooFit2(TTree *tree) int level = 3; std::cout << "num entries = " << data.numEntries() << std::endl; bool save = true; - (pdf[N - 1]->getVariables())->Print("v"); // print the parameters + std::unique_ptr{pdf[N - 1]->getVariables()}->Print("v"); // print the parameters std::cout << "\n\nDo the fit now \n\n"; #else int level = -1; diff --git a/roofit/roofit/test/vectorisedPDFs/VectorisedPDFTests.cxx b/roofit/roofit/test/vectorisedPDFs/VectorisedPDFTests.cxx index 6d2c5adb01e7d..282865e692304 100644 --- a/roofit/roofit/test/vectorisedPDFs/VectorisedPDFTests.cxx +++ b/roofit/roofit/test/vectorisedPDFs/VectorisedPDFTests.cxx @@ -457,7 +457,7 @@ std::unique_ptr PDFTest::runBatchFit(RooAbsPdf *pdf) kickParameters(); makePlots(::testing::UnitTest::GetInstance()->current_test_info()->name() + std::string("_batch_prefit")); - auto pars = pdf->getParameters(*_dataFit); + std::unique_ptr pars{pdf->getParameters(*_dataFit)}; *pars = _parameters; for (unsigned int index = 0; index < pars->size(); ++index) { diff --git a/roofit/roofitcore/src/RooAbsReal.cxx b/roofit/roofitcore/src/RooAbsReal.cxx index da039f60240cd..0c238a1861b12 100644 --- a/roofit/roofitcore/src/RooAbsReal.cxx +++ b/roofit/roofitcore/src/RooAbsReal.cxx @@ -2247,10 +2247,10 @@ RooPlot* RooAbsReal::plotAsymOn(RooPlot *frame, const RooAbsCategoryLValue& asym // Take out data-projected dependents from projectedVars - RooArgSet* projDataNeededVars = nullptr ; + std::unique_ptr projDataNeededVars; if (o.projData) { - projDataNeededVars = projectedVars.selectCommon(projDataVars); - projectedVars.remove(projDataVars,true,true) ; + projDataNeededVars.reset(projectedVars.selectCommon(projDataVars)); + projectedVars.remove(projDataVars, true, true); } // Take out plotted asymmetry from projection diff --git a/roofit/roofitcore/src/RooRealMPFE.cxx b/roofit/roofitcore/src/RooRealMPFE.cxx index cd87846508e01..cc3529a15adb9 100644 --- a/roofit/roofitcore/src/RooRealMPFE.cxx +++ b/roofit/roofitcore/src/RooRealMPFE.cxx @@ -53,6 +53,7 @@ For general multiprocessing in ROOT, please refer to the TProcessExecutor class. #endif #include +#include #include #include "RooRealMPFE.h" #include "RooArgSet.h" @@ -170,7 +171,7 @@ void RooRealMPFE::initVars() _saveVars.removeAll() ; // Retrieve non-constant parameters - auto vars = _arg->getParameters(RooArgSet()); + std::unique_ptr vars{_arg->getParameters(RooArgSet())}; // RooArgSet *ncVars = vars->selectByAttrib("Constant", false); RooArgList varList(*vars) ; diff --git a/roofit/roofitcore/src/RooVectorDataStore.cxx b/roofit/roofitcore/src/RooVectorDataStore.cxx index 96c20887e309c..2f1b22d1b64b8 100644 --- a/roofit/roofitcore/src/RooVectorDataStore.cxx +++ b/roofit/roofitcore/src/RooVectorDataStore.cxx @@ -829,6 +829,7 @@ void RooVectorDataStore::cacheArgs(const RooAbsArg* owner, RooArgSet& newVarSet, std::vector nsetList ; std::vector> argObsList ; + std::vector> ownedNsets; // Now need to attach branch buffers of clones for (const auto arg : cloneSet) { @@ -844,8 +845,8 @@ void RooVectorDataStore::cacheArgs(const RooAbsArg* owner, RooArgSet& newVarSet, // std::cout << "RooVectorDataStore::cacheArgs() cached node " << arg->GetName() << " has a normalization set specification CATNormSet = " << catNset << std::endl ; RooArgSet anset = RooHelpers::selectFromArgSet(nset ? *nset : RooArgSet{}, catNset); - normSet = anset.selectCommon(*argObs); - + ownedNsets.emplace_back(anset.selectCommon(*argObs)); + normSet = ownedNsets.back().get(); } const char* catCset = arg->getStringAttribute("CATCondSet") ; if (catCset) { diff --git a/roofit/roofitcore/src/TestStatistics/RooAbsL.cxx b/roofit/roofitcore/src/TestStatistics/RooAbsL.cxx index 76ad278c4dfc8..badd019b14964 100644 --- a/roofit/roofitcore/src/TestStatistics/RooAbsL.cxx +++ b/roofit/roofitcore/src/TestStatistics/RooAbsL.cxx @@ -129,7 +129,7 @@ void RooAbsL::initClones(RooAbsPdf &inpdf, RooAbsData &indata) // ****************************************************************** // Attach FUNC to data set - auto _funcObsSet = pdf_->getObservables(indata); + std::unique_ptr _funcObsSet{pdf_->getObservables(indata)}; if (pdf_->getAttribute("BinnedLikelihood")) { pdf_->setAttribute("BinnedLikelihoodActive"); diff --git a/roofit/roostats/src/AsymptoticCalculator.cxx b/roofit/roostats/src/AsymptoticCalculator.cxx index 450604729bee8..899c604faeb1d 100644 --- a/roofit/roostats/src/AsymptoticCalculator.cxx +++ b/roofit/roostats/src/AsymptoticCalculator.cxx @@ -245,8 +245,8 @@ bool AsymptoticCalculator::Initialize() const { if (!fNominalAsimov) { if (verbose >= 0) oocoutI(nullptr,InputArguments) << "AsymptoticCalculator: Asimov data will be generated using fitted nuisance parameter values" << std::endl; - RooArgSet * tmp = (RooArgSet*) poiAlt.snapshot(); - fAsimovData = MakeAsimovData( data, *GetNullModel(), poiAlt, fAsimovGlobObs,tmp); + std::unique_ptr tmp{static_cast(poiAlt.snapshot())}; + fAsimovData = MakeAsimovData(data, *GetNullModel(), poiAlt, fAsimovGlobObs, tmp.get()); } else { diff --git a/roofit/roostats/src/HypoTestInverter.cxx b/roofit/roostats/src/HypoTestInverter.cxx index 9fa2bed777aaf..16f12bb86d21a 100644 --- a/roofit/roostats/src/HypoTestInverter.cxx +++ b/roofit/roostats/src/HypoTestInverter.cxx @@ -1237,7 +1237,8 @@ SamplingDistribution * HypoTestInverter::RebuildDistributions(bool isUpper, int inverter.SetData(*bkgdata); // print global observables - auto gobs = bModel->GetPdf()->getVariables()->selectCommon(* sbModel->GetGlobalObservables() ); + std::unique_ptr bModelVars{bModel->GetPdf()->getVariables()}; + std::unique_ptr gobs{bModelVars->selectCommon(*sbModel->GetGlobalObservables())}; gobs->Print("v"); HypoTestInverterResult * r = inverter.GetInterval(); diff --git a/roofit/roostats/src/LikelihoodIntervalPlot.cxx b/roofit/roostats/src/LikelihoodIntervalPlot.cxx index 5d4d669f5aa28..d072e35d1459a 100644 --- a/roofit/roostats/src/LikelihoodIntervalPlot.cxx +++ b/roofit/roostats/src/LikelihoodIntervalPlot.cxx @@ -194,7 +194,8 @@ void LikelihoodIntervalPlot::Draw(const Option_t *options) const double xcont_min = fInterval->LowerLimit(*myparam); const double xcont_max = fInterval->UpperLimit(*myparam); - RooRealVar* myarg = static_cast(newProfile->getVariables()->find(myparam->GetName())); + std::unique_ptr vars{newProfile->getVariables()}; + RooRealVar *myarg = static_cast(vars->find(myparam->GetName())); double x1 = myarg->getMin(); double x2 = myarg->getMax(); @@ -336,7 +337,8 @@ void LikelihoodIntervalPlot::Draw(const Option_t *options) double cont_level = ROOT::Math::chisquared_quantile(fInterval->ConfidenceLevel(),fNdimPlot); // level for -2log LR cont_level = cont_level/2; // since we are plotting -log LR - RooArgList params(*newProfile->getVariables()); + std::unique_ptr vars{newProfile->getVariables()}; + RooArgList params(*vars); // set values and error for the POI to the best fit values for (std::size_t i = 0; i < params.size(); ++i) { RooRealVar & par = static_cast( params[i]); diff --git a/roofit/roostats/src/PdfProposal.cxx b/roofit/roostats/src/PdfProposal.cxx index bbfbd2d0199f9..7abba52c9235e 100644 --- a/roofit/roostats/src/PdfProposal.cxx +++ b/roofit/roostats/src/PdfProposal.cxx @@ -177,8 +177,9 @@ double PdfProposal::GetProposalDensity(RooArgSet& x1, RooArgSet& x2) void PdfProposal::AddMapping(RooRealVar& proposalParam, RooAbsReal& update) { - fMaster.add(*update.getParameters(static_cast(nullptr))); - if (update.getParameters(static_cast(nullptr))->empty()) + std::unique_ptr params{update.getParameters(static_cast(nullptr))}; + fMaster.add(*params); + if (params->empty()) fMaster.add(update); fMap.insert(std::pair(&proposalParam, &update)); } diff --git a/roofit/roostats/src/SPlot.cxx b/roofit/roostats/src/SPlot.cxx index c18eefd053ecd..227948497bcf0 100644 --- a/roofit/roostats/src/SPlot.cxx +++ b/roofit/roostats/src/SPlot.cxx @@ -466,7 +466,8 @@ void SPlot::AddSWeight(RooAbsPdf *pdf, const RooArgList &yieldsTmp, const RooArg } const Int_t nspec = yieldsTmp.size(); - RooArgList yields = *static_cast(yieldsTmp.snapshot(false)); + std::unique_ptr yieldsSnapshot{static_cast(yieldsTmp.snapshot(false))}; + RooArgList &yields = *yieldsSnapshot; if (RooMsgService::instance().isActive(this, RooFit::InputArguments, RooFit::DEBUG)) { coutI(InputArguments) << "Printing Yields" << std::endl; diff --git a/roofit/roostats/src/SimpleLikelihoodRatioTestStat.cxx b/roofit/roostats/src/SimpleLikelihoodRatioTestStat.cxx index 03731ee4a1135..324514445e2c3 100644 --- a/roofit/roostats/src/SimpleLikelihoodRatioTestStat.cxx +++ b/roofit/roostats/src/SimpleLikelihoodRatioTestStat.cxx @@ -37,8 +37,10 @@ double RooStats::SimpleLikelihoodRatioTestStat::Evaluate(RooAbsData& data, RooAr // strip pdfs of constraints (which cancel out in the ratio) to avoid unnecessary computations and computational errors if (fFirstEval) { - fNullPdf = RooStats::MakeUnconstrainedPdf(*fNullPdf, *fNullPdf->getObservables(data)); - fAltPdf = RooStats::MakeUnconstrainedPdf(*fAltPdf , *fAltPdf->getObservables(data) ); + std::unique_ptr nullObs{fNullPdf->getObservables(data)}; + fNullPdf = RooStats::MakeUnconstrainedPdf(*fNullPdf, *nullObs); + std::unique_ptr altObs{fAltPdf->getObservables(data)}; + fAltPdf = RooStats::MakeUnconstrainedPdf(*fAltPdf, *altObs); } fFirstEval = false; diff --git a/roofit/roostats/src/ToyMCImportanceSampler.cxx b/roofit/roostats/src/ToyMCImportanceSampler.cxx index 2c997814b52f5..377c3a09699db 100644 --- a/roofit/roostats/src/ToyMCImportanceSampler.cxx +++ b/roofit/roostats/src/ToyMCImportanceSampler.cxx @@ -302,7 +302,7 @@ RooAbsData* ToyMCImportanceSampler::GenerateToyData( std::unique_ptr allVarsImpDens{fImportanceDensities[fIndexGenDensity]->getVariables()}; allVars->add(*allVarsImpDens); } - const RooArgSet* saveVars = (const RooArgSet*)allVars->snapshot(); + std::unique_ptr saveVars{static_cast(allVars->snapshot())}; double globalWeight = 1.0; if(fNuisanceParametersSampler) { // use nuisance parameters? @@ -413,10 +413,7 @@ RooAbsData* ToyMCImportanceSampler::GenerateToyData( ooccoutD(nullptr,InputArguments) << "weights["<assign(*saveVars); - delete saveVars; return data; } diff --git a/roofit/xroofit/src/xRooNode.cxx b/roofit/xroofit/src/xRooNode.cxx index 72059cac8c29e..bd856470a5a5d 100644 --- a/roofit/xroofit/src/xRooNode.cxx +++ b/roofit/xroofit/src/xRooNode.cxx @@ -9016,9 +9016,9 @@ TH1 *xRooNode::BuildHistogram(RooAbsLValue *v, bool empty, bool errors, int binS for (auto pdf : bins()) { // auto _pdf = // pdf->get()->createProjection(*pdf->get()->getObservables(*_obs.get())); - auto _pdf = - new xRooProjectedPdf(TString::Format("%s_projection", pdf->GetName()), "", *pdf->get(), - *pdf->get()->getObservables(*_obs.get())); + std::unique_ptr projObs{pdf->get()->getObservables(*_obs.get())}; + auto _pdf = new xRooProjectedPdf(TString::Format("%s_projection", pdf->GetName()), "", + *pdf->get(), *projObs); if (hasRange) { dynamic_cast(_pdf)->setNormRange("coordRange"); } From 83d5f970337c37ee0c298cfe6e3e7a27c2feb94c Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Fri, 19 Jun 2026 08:13:11 +0200 Subject: [PATCH 3/5] [RF] Modernize RooFit code with static_range_cast and dynamic_range_cast The `static_range_cast` and `dynamic_range_cast` are RooFit helper functions to iterate over polymorphic collections like `RooArgSet` and `RooArgList`. Using them instead of casting in a separate line makes the code more concise with fewer lines and throw-away temoprary variables. (cherry picked from commit 4f2cf43e9fa2bd576087a17c9d5a6ccffd7358db) --- roofit/hs3/src/JSONFactories_HistFactory.cxx | 3 +- roofit/hs3/src/JSONFactories_RooFitCore.cxx | 3 +- roofit/roofit/src/RooJeffreysPrior.cxx | 7 +-- roofit/roofit/src/RooLagrangianMorphFunc.cxx | 62 ++++++------------- roofit/roofit/src/RooLegacyExpPoly.cxx | 7 +-- roofit/roofit/src/RooMultiBinomial.cxx | 4 +- roofit/roofit/src/RooNDKeysPdf.cxx | 8 +-- roofit/roofitcore/src/FitHelpers.cxx | 7 +-- roofit/roofitcore/src/RooAbsAnaConvPdf.cxx | 6 +- roofit/roofitcore/src/RooAddModel.cxx | 27 +++----- roofit/roofitcore/src/RooAddPdf.cxx | 8 +-- roofit/roofitcore/src/RooClassFactory.cxx | 4 +- roofit/roofitcore/src/RooFFTConvPdf.cxx | 3 +- roofit/roofitcore/src/RooFoamGenerator.cxx | 3 +- .../roofitcore/src/RooLinearCombination.cxx | 8 +-- roofit/roofitcore/src/RooSuperCategory.cxx | 3 +- roofit/roostats/src/BayesianCalculator.cxx | 3 +- .../roostats/src/DetailedOutputAggregator.cxx | 4 +- roofit/roostats/src/HypoTestInverter.cxx | 3 +- roofit/roostats/src/LikelihoodInterval.cxx | 9 ++- .../roostats/src/LikelihoodIntervalPlot.cxx | 7 +-- .../src/ProfileLikelihoodCalculator.cxx | 9 ++- 22 files changed, 76 insertions(+), 122 deletions(-) diff --git a/roofit/hs3/src/JSONFactories_HistFactory.cxx b/roofit/hs3/src/JSONFactories_HistFactory.cxx index df59000016fc0..f3cffdbb7a9f5 100644 --- a/roofit/hs3/src/JSONFactories_HistFactory.cxx +++ b/roofit/hs3/src/JSONFactories_HistFactory.cxx @@ -1412,8 +1412,7 @@ class HistFactoryStreamer_ProdPdf : public RooFit::JSONIO::Exporter { { std::vector constraints; RooRealSumPdf *sumpdf = nullptr; - for (RooAbsArg *v : prodpdf->pdfList()) { - RooAbsPdf *pdf = static_cast(v); + for (auto *pdf : static_range_cast(prodpdf->pdfList())) { auto thispdf = dynamic_cast(pdf); if (thispdf) { if (!sumpdf) diff --git a/roofit/hs3/src/JSONFactories_RooFitCore.cxx b/roofit/hs3/src/JSONFactories_RooFitCore.cxx index 90ffd950c7fe6..29a3f66fa99d4 100644 --- a/roofit/hs3/src/JSONFactories_RooFitCore.cxx +++ b/roofit/hs3/src/JSONFactories_RooFitCore.cxx @@ -627,8 +627,7 @@ class ParamHistFuncFactory : public RooFit::JSONIO::Importer { // Now build the final list following the order in varList RooArgList vars; - for (int i = 0; i < varList.getSize(); ++i) { - const auto *refVar = dynamic_cast(varList.at(i)); + for (auto *refVar : dynamic_range_cast(varList)) { if (!refVar) continue; diff --git a/roofit/roofit/src/RooJeffreysPrior.cxx b/roofit/roofit/src/RooJeffreysPrior.cxx index 50d1e451c140d..bb907e4046d0d 100644 --- a/roofit/roofit/src/RooJeffreysPrior.cxx +++ b/roofit/roofit/src/RooJeffreysPrior.cxx @@ -82,11 +82,10 @@ double RooJeffreysPrior::evaluate() const auto& pdf = _nominal.arg(); RooAbsPdf* clonePdf = static_cast(pdf.cloneTree()); std::unique_ptr vars{clonePdf->getParameters(_obsSet)}; - for (auto varTmp : *vars) { - auto& var = static_cast(*varTmp); - auto range = var.getRange(); + for (auto* var : static_range_cast(*vars)) { + auto range = var->getRange(); double span = range.second - range.first; - var.setRange(range.first - 0.1*span, range.second + 0.1 * span); + var->setRange(range.first - 0.1*span, range.second + 0.1 * span); } cacheElm = new CacheElem; diff --git a/roofit/roofit/src/RooLagrangianMorphFunc.cxx b/roofit/roofit/src/RooLagrangianMorphFunc.cxx index 94473baec0a1f..887eba52ea8b1 100644 --- a/roofit/roofit/src/RooLagrangianMorphFunc.cxx +++ b/roofit/roofit/src/RooLagrangianMorphFunc.cxx @@ -461,8 +461,7 @@ void setOwnerRecursive(TFolder *theFolder) theFolder->SetOwner(); // And also need to set up ownership for nested folders auto subdirs = theFolder->GetListOfFolders(); - for (auto *subdir : *subdirs) { - auto thisfolder = dynamic_cast(subdir); + for (auto *thisfolder : dynamic_range_cast(*subdirs)) { if (thisfolder) { // no explicit deletion here, will be handled by parent setOwnerRecursive(thisfolder); @@ -678,8 +677,7 @@ inline bool setParam(RooRealVar *p, double val, bool force) template inline bool setParams(const T2 &args, T1 val) { - for (auto itr : args) { - RooRealVar *param = dynamic_cast(itr); + for (auto *param : dynamic_range_cast(args)) { if (!param) continue; setParam(param, val, true); @@ -696,8 +694,7 @@ inline bool setParams(const std::map &point, const T2 &args, bool force = false, T1 defaultVal = 0) { bool ok = true; - for (auto itr : args) { - RooRealVar *param = dynamic_cast(itr); + for (auto *param : dynamic_range_cast(args)) { if (!param || param->isConstant()) continue; ok = setParam(param, defaultVal, force) && ok; @@ -725,8 +722,7 @@ inline bool setParams(TH1 *hist, const T &args, bool force = false) { bool ok = true; - for (auto itr : args) { - RooRealVar *param = dynamic_cast(itr); + for (auto *param : dynamic_range_cast(args)) { if (!param) continue; ok = setParam(param, 0., force) && ok; @@ -752,8 +748,7 @@ template inline RooLagrangianMorphFunc::ParamSet getParams(const T ¶meters) { RooLagrangianMorphFunc::ParamSet retval; - for (auto itr : parameters) { - RooRealVar *param = dynamic_cast(itr); + for (auto *param : dynamic_range_cast(parameters)) { if (!param) continue; retval[param->GetName()] = param->getVal(); @@ -977,9 +972,7 @@ inline void fillFeynmanDiagram(FeynmanDiagram &diagram, const std::vector vertexCouplings(ncouplings, false); int idx = -1; - RooAbsReal *coupling; - for (auto citr : couplings) { - coupling = dynamic_cast(citr); + for (auto *coupling : dynamic_range_cast(couplings)) { idx++; if (!coupling) { std::cerr << "encountered invalid list of couplings in vertex!" << std::endl; @@ -1120,11 +1113,9 @@ FormulaList buildFormulas(const char *mfname, const RooLagrangianMorphFunc::Para std::cerr << "internal error, number of operators inconsistent!" << std::endl; } - RooAbsReal *obj0; int idx = 0; - for (auto itr1 : couplings) { - obj0 = dynamic_cast(itr1); + for (auto *obj0 : dynamic_range_cast(couplings)) { if (obj0->getVal() != 0) { couplingsZero[idx] = false; } @@ -1132,8 +1123,7 @@ FormulaList buildFormulas(const char *mfname, const RooLagrangianMorphFunc::Para } } - for (auto itr2 : flags) { - auto obj1 = dynamic_cast(itr2); + for (auto *obj1 : dynamic_range_cast(flags)) { int nZero = 0; int nNonZero = 0; for (auto sampleit : inputFlags) { @@ -1203,8 +1193,7 @@ FormulaList buildFormulas(const char *mfname, const RooLagrangianMorphFunc::Para // check and apply flags bool removedByFlag = false; - for (auto itr : flags) { - auto obj = dynamic_cast(itr); + for (auto *obj : dynamic_range_cast(flags)) { if (!obj) continue; TString sval(obj->getStringAttribute("NewPhysics")); @@ -1411,10 +1400,8 @@ class RooLagrangianMorphFunc::CacheElem : public RooAbsCacheElement { // set all vars to value stored in input file setParams(sampleit.second, operators, true); bool first = true; - RooAbsReal *obj; - for (auto itr : _couplings) { - obj = dynamic_cast(itr); + for (auto *obj : dynamic_range_cast(_couplings)) { if (!first) std::cerr << ", "; oocxcoutW((TObject *)nullptr, Eval) << obj->GetName() << "=" << obj->getVal(); @@ -1953,8 +1940,8 @@ RooLagrangianMorphFunc::RooLagrangianMorphFunc(const RooLagrangianMorphFunc &oth { for (size_t j = 0; j < other._diagrams.size(); ++j) { std::vector diagram; - for (size_t i = 0; i < other._diagrams[j].size(); ++i) { - RooListProxy *list = new RooListProxy(other._diagrams[j][i]->GetName(), this, *(other._diagrams[j][i])); + for (auto *elem : other._diagrams[j]) { + RooListProxy *list = new RooListProxy(elem->GetName(), this, *elem); diagram.push_back(list); } _diagrams.push_back(diagram); @@ -2158,8 +2145,7 @@ RooProduct *RooLagrangianMorphFunc::getSumElement(const char *name) const prodname.Append("_"); prodname.Append(this->GetName()); - for (auto itr : *args) { - RooProduct *prod = dynamic_cast(itr); + for (auto *prod : dynamic_range_cast(*args)) { if (!prod) continue; TString sname(prod->GetName()); @@ -2215,11 +2201,9 @@ void RooLagrangianMorphFunc::printSampleWeights() const void RooLagrangianMorphFunc::randomizeParameters(double z) { - RooRealVar *obj; TRandom3 r; - for (auto itr : _operators) { - obj = dynamic_cast(itr); + for (auto *obj : dynamic_range_cast(_operators)) { double val = obj->getVal(); if (obj->isConstant()) continue; @@ -2504,8 +2488,7 @@ RooLagrangianMorphFunc::ParamSet RooLagrangianMorphFunc::getMorphParameters(cons void RooLagrangianMorphFunc::setParameters(const RooArgList *list) { - for (auto itr : *list) { - RooRealVar *param = dynamic_cast(itr); + for (auto *param : dynamic_range_cast(*list)) { if (!param) continue; this->setParameter(param->GetName(), param->getVal()); @@ -2562,8 +2545,7 @@ TH1 *RooLagrangianMorphFunc::createTH1(const std::string &name, bool correlateEr double val = 0; double unc2 = 0; double unc = 0; - for (auto itr : *args) { - RooProduct *prod = dynamic_cast(itr); + for (auto *prod : dynamic_range_cast(*args)) { if (!prod) continue; RooAbsArg *phys = prod->components().find(Form("phys_%s", prod->GetName())); @@ -2595,8 +2577,7 @@ int RooLagrangianMorphFunc::countContributingFormulas() const if (!mf) coutE(InputArguments) << "unable to retrieve morphing function" << std::endl; std::unique_ptr args{mf->getComponents()}; - for (auto itr : *args) { - RooProduct *prod = dynamic_cast(itr); + for (auto *prod : dynamic_range_cast(*args)) { if (prod->getVal() != 0) { nFormulas++; } @@ -2721,8 +2702,7 @@ const RooArgList *RooLagrangianMorphFunc::getCouplingSet() const RooLagrangianMorphFunc::ParamSet RooLagrangianMorphFunc::getCouplings() const { RooLagrangianMorphFunc::ParamSet couplings; - for (auto obj : *(this->getCouplingSet())) { - RooAbsReal *var = dynamic_cast(obj); + for (auto *var : dynamic_range_cast(*(this->getCouplingSet()))) { if (!var) continue; const std::string name(var->GetName()); @@ -2847,8 +2827,7 @@ double RooLagrangianMorphFunc::expectedUncertainty() const void RooLagrangianMorphFunc::printParameters() const { // print the parameters and their current values - for (auto obj : _operators) { - RooRealVar *param = static_cast(obj); + for (auto *param : static_range_cast(_operators)) { if (!param) continue; param->Print(); @@ -2860,8 +2839,7 @@ void RooLagrangianMorphFunc::printParameters() const void RooLagrangianMorphFunc::printFlags() const { - for (auto flag : _flags) { - RooRealVar *param = static_cast(flag); + for (auto *param : static_range_cast(_flags)) { if (!param) continue; param->Print(); diff --git a/roofit/roofit/src/RooLegacyExpPoly.cxx b/roofit/roofit/src/RooLegacyExpPoly.cxx index 8a28861393f41..af79382b7a1f1 100644 --- a/roofit/roofit/src/RooLegacyExpPoly.cxx +++ b/roofit/roofit/src/RooLegacyExpPoly.cxx @@ -109,8 +109,8 @@ double RooLegacyExpPoly::evaluateLog() const const double x = _x; double xpow = std::pow(x, lowestOrder); double retval = 0; - for (size_t i = 0; i < sz; ++i) { - retval += coefs[i] * xpow; + for (double coef : coefs) { + retval += coef * xpow; xpow *= x; } @@ -158,9 +158,8 @@ void RooLegacyExpPoly::adjustLimits() if (x) { const double xmax = x->getMax(); double xmaxpow = std::pow(xmax, lowestOrder); - for (size_t i = 0; i < sz; ++i) { + for (auto *coef : dynamic_range_cast(_coefList)) { double thismax = max / xmaxpow; - RooRealVar *coef = dynamic_cast(this->_coefList.at(i)); if (coef) { coef->setVal(thismax); coef->setMax(thismax); diff --git a/roofit/roofit/src/RooMultiBinomial.cxx b/roofit/roofit/src/RooMultiBinomial.cxx index 2699a7afdb310..adcff7d05cbed 100644 --- a/roofit/roofit/src/RooMultiBinomial.cxx +++ b/roofit/roofit/src/RooMultiBinomial.cxx @@ -126,8 +126,8 @@ double RooMultiBinomial::evaluate() const // Calculate efficiency for combination of accept/reject categories // put equal to zero if combination of only zeros AND chosen to be invisible - for (int i=0; inEventsBMSW; if (norm<0.) norm=0.; - for (Int_t i=0; isIdcs.size()); ++i) { + for (Int_t sIdx : bi->sIdcs) { double prob=1.; - const vector& x = _dataPts[bi->sIdcs[i]]; - const vector& weight = (*_weights)[_idx[bi->sIdcs[i]]]; + const vector& x = _dataPts[sIdx]; + const vector& weight = (*_weights)[_idx[sIdx]]; vector chi(_nDim,100.); @@ -1202,7 +1202,7 @@ double RooNDKeysPdf::analyticalIntegral(Int_t code, const char* rangeName) const } } - norm += prob * _wMap.at(_idx[bi->sIdcs[i]]); + norm += prob * _wMap.at(_idx[sIdx]); } cxcoutD(Eval) << "RooNDKeysPdf::analyticalIntegral() : Final normalization : " << norm << " " << bi->nEventsBW << std::endl; diff --git a/roofit/roofitcore/src/FitHelpers.cxx b/roofit/roofitcore/src/FitHelpers.cxx index e1f1b3512a11e..0c53e2073f50a 100644 --- a/roofit/roofitcore/src/FitHelpers.cxx +++ b/roofit/roofitcore/src/FitHelpers.cxx @@ -784,8 +784,7 @@ std::unique_ptr createNLL(RooAbsPdf &pdf, RooAbsData &data, const Ro // Create range with name 'fit' with above limits on all observables RooArgSet obs; pdf.getObservables(data.get(), obs); - for (auto arg : obs) { - RooRealVar *rrv = dynamic_cast(arg); + for (auto *rrv : dynamic_range_cast(obs)) { if (rrv) rrv->setRange("fit", rangeLo, rangeHi); } @@ -1043,8 +1042,8 @@ std::unique_ptr createChi2FromDataHist(RooAbsReal &real, RooDataHist const double rangeHi = pc.getDouble("rangeHi"); RooArgSet obs; real.getObservables(data.get(), obs); - for (auto arg : obs) { - if (auto *rrv = dynamic_cast(arg)) { + for (auto *rrv : dynamic_range_cast(obs)) { + if (rrv) { rrv->setRange("fit", rangeLo, rangeHi); } } diff --git a/roofit/roofitcore/src/RooAbsAnaConvPdf.cxx b/roofit/roofitcore/src/RooAbsAnaConvPdf.cxx index 85fe698dab089..d9c9dffd50277 100644 --- a/roofit/roofitcore/src/RooAbsAnaConvPdf.cxx +++ b/roofit/roofitcore/src/RooAbsAnaConvPdf.cxx @@ -209,8 +209,7 @@ bool RooAbsAnaConvPdf::changeModel(const RooResolutionModel& newModel) { RooArgList newConvSet ; bool allOK(true) ; - for (auto convArg : _convSet) { - auto conv = static_cast(convArg); + for (auto *conv : static_range_cast(_convSet)) { // Build new resolution model std::unique_ptr newConv{newModel.convolution(const_cast(&conv->basis()),this)}; @@ -330,8 +329,7 @@ double RooAbsAnaConvPdf::evaluate() const double result(0) ; Int_t index(0) ; - for (auto convArg : _convSet) { - auto conv = static_cast(convArg); + for (auto *conv : static_range_cast(_convSet)) { double coef = coefficient(index++) ; if (coef!=0.) { const double c = conv->getVal(nullptr); diff --git a/roofit/roofitcore/src/RooAddModel.cxx b/roofit/roofitcore/src/RooAddModel.cxx index 2a33a9ea3a1c2..ba7b3c181c9c3 100644 --- a/roofit/roofitcore/src/RooAddModel.cxx +++ b/roofit/roofitcore/src/RooAddModel.cxx @@ -244,8 +244,7 @@ RooResolutionModel* RooAddModel::convolution(RooFormulaVar* inBasis, RooAbsArg* newTitle.Append(inBasis->GetName()) ; RooArgList modelList ; - for (auto obj : _pdfList) { - auto model = static_cast(obj); + for (auto *model : static_range_cast(_pdfList)) { // Create component convolution RooResolutionModel* conv = model->convolution(inBasis,owner) ; modelList.add(*conv) ; @@ -281,8 +280,7 @@ Int_t RooAddModel::basisCode(const char* name) const { bool first(true); bool code(false); - for (auto obj : _pdfList) { - auto model = static_cast(obj); + for (auto *model : static_range_cast(_pdfList)) { Int_t subCode = model->basisCode(name) ; if (first) { code = subCode ; @@ -362,8 +360,7 @@ double RooAddModel::evaluate() const double snormVal ; double value(0) ; Int_t i(0) ; - for (auto obj : _pdfList) { - auto pdf = static_cast(obj); + for (auto *pdf : static_range_cast(_pdfList)) { if (_coefCache[i]!=0.) { snormVal = nset ? cache->suppNormVal(i) : 1.0 ; @@ -496,8 +493,7 @@ void RooAddModel::getCompIntList(const RooArgSet* nset, const RooArgSet* iset, p cache = new IntCacheElem ; // Fill Cache - for (auto obj : _pdfList) { - auto model = static_cast(obj); + for (auto *model : static_range_cast(_pdfList)) { cache->_intList.addOwned(std::unique_ptr{model->createIntegral(*iset,nset,nullptr,isetRangeName)}); } @@ -550,8 +546,7 @@ double RooAddModel::analyticalIntegralWN(Int_t code, const RooArgSet* normSet, c double snormVal ; double value(0) ; Int_t i(0) ; - for (const auto obj : *compIntList) { - auto pdfInt = static_cast(obj); + for (auto *pdfInt : static_range_cast(*compIntList)) { if (_coefCache[i]!=0.) { snormVal = nset ? pcache->suppNormVal(i) : 1.0 ; double intVal = pdfInt->getVal(nset) ; @@ -579,16 +574,14 @@ double RooAddModel::expectedEvents(const RooArgSet* nset) const if (_allExtendable) { // Sum of the extended terms - for (auto obj : _pdfList) { - auto pdf = static_cast(obj); + for (auto *pdf : static_range_cast(_pdfList)) { expectedTotal += pdf->expectedEvents(nset) ; } } else { // Sum the coefficients - for (const auto obj : _coefList) { - auto coef = static_cast(obj); + for (auto *coef : static_range_cast(_coefList)) { expectedTotal += coef->getVal() ; } } @@ -651,8 +644,7 @@ RooAbsGenContext* RooAddModel::genContext(const RooArgSet &vars, const RooDataSe bool RooAddModel::isDirectGenSafe(const RooAbsArg& arg) const { - for (auto obj : _pdfList) { - auto pdf = static_cast(obj); + for (auto *pdf : static_range_cast(_pdfList)) { if (!pdf->isDirectGenSafe(arg)) { return false ; @@ -668,8 +660,7 @@ bool RooAddModel::isDirectGenSafe(const RooAbsArg& arg) const Int_t RooAddModel::getGenerator(const RooArgSet& directVars, RooArgSet &/*generateVars*/, bool /*staticInitOK*/) const { - for (auto obj : _pdfList) { - auto pdf = static_cast(obj); + for (auto *pdf : static_range_cast(_pdfList)) { RooArgSet tmp ; if (pdf->getGenerator(directVars,tmp)==0) { diff --git a/roofit/roofitcore/src/RooAddPdf.cxx b/roofit/roofitcore/src/RooAddPdf.cxx index edc8f8963ec79..d55e89e0f3aeb 100644 --- a/roofit/roofitcore/src/RooAddPdf.cxx +++ b/roofit/roofitcore/src/RooAddPdf.cxx @@ -774,12 +774,12 @@ double RooAddPdf::expectedEvents(const RooArgSet* nset) const } else { if (_allExtendable) { - for(auto const& arg : _pdfList) { - expectedTotal += static_cast(arg)->expectedEvents(nset) ; + for(auto *arg : static_range_cast(_pdfList)) { + expectedTotal += arg->expectedEvents(nset) ; } } else { - for(auto const& arg : _coefList) { - expectedTotal += static_cast(arg)->getVal(nset) ; + for(auto *arg : static_range_cast(_coefList)) { + expectedTotal += arg->getVal(nset) ; } } diff --git a/roofit/roofitcore/src/RooClassFactory.cxx b/roofit/roofitcore/src/RooClassFactory.cxx index ffa0e13920aa8..65e4179edc037 100644 --- a/roofit/roofitcore/src/RooClassFactory.cxx +++ b/roofit/roofitcore/src/RooClassFactory.cxx @@ -382,9 +382,9 @@ std::string listVars(std::vector const &alist, std::vector co std::string declareVarSpans(std::vector const &alist) { std::stringstream ss; - for (std::size_t i = 0; i < alist.size(); ++i) { + for (auto const &elem : alist) { ss << " " - << "std::span " << alist[i] << "Span = ctx.at(" << alist[i] << ");\n"; + << "std::span " << elem << "Span = ctx.at(" << elem << ");\n"; } return ss.str(); } diff --git a/roofit/roofitcore/src/RooFFTConvPdf.cxx b/roofit/roofitcore/src/RooFFTConvPdf.cxx index 790ea9aef1f5b..c6895b6d4936f 100644 --- a/roofit/roofitcore/src/RooFFTConvPdf.cxx +++ b/roofit/roofitcore/src/RooFFTConvPdf.cxx @@ -520,8 +520,7 @@ void RooFFTConvPdf::fillCacheObject(RooAbsCachedPdf::PdfCacheElem& cache) const std::vector obsLV(n); Int_t i(0) ; - for (auto const& arg : otherObs) { - RooAbsLValue* lvarg = dynamic_cast(arg) ; + for (auto* lvarg : dynamic_range_cast(otherObs)) { obsLV[i] = lvarg ; binCur[i] = 0 ; // coverity[FORWARD_NULL] diff --git a/roofit/roofitcore/src/RooFoamGenerator.cxx b/roofit/roofitcore/src/RooFoamGenerator.cxx index 11369d7c86bd4..98709e66845f6 100644 --- a/roofit/roofitcore/src/RooFoamGenerator.cxx +++ b/roofit/roofitcore/src/RooFoamGenerator.cxx @@ -149,8 +149,7 @@ const RooArgSet *RooFoamGenerator::generateEvent(UInt_t /*remaining*/, double& / // Transfer contents to dataset Int_t i(0) ; - for (auto arg : _realVars) { - auto var = static_cast(arg); + for (auto* var : static_range_cast(_realVars)) { var->setVal(_xmin[i] + _range[i]*_vec[i]) ; i++ ; } diff --git a/roofit/roofitcore/src/RooLinearCombination.cxx b/roofit/roofitcore/src/RooLinearCombination.cxx index e030befd5bb5c..fd0b95d9c02cf 100644 --- a/roofit/roofitcore/src/RooLinearCombination.cxx +++ b/roofit/roofitcore/src/RooLinearCombination.cxx @@ -133,8 +133,8 @@ std::list *RooLinearCombination::binBoundaries(RooAbsRealLValue &obs, double xhi) const { // Forward the plot sampling hint from the p.d.f. that defines the observable // obs - for(auto const& func : _actualVars) { - auto binb = static_cast(func)->binBoundaries(obs, xlo, xhi); + for(auto *func : static_range_cast(_actualVars)) { + auto binb = func->binBoundaries(obs, xlo, xhi); if (binb) { return binb; } @@ -147,8 +147,8 @@ std::list *RooLinearCombination::plotSamplingHint(RooAbsRealLValue &obs, double xhi) const { // Forward the plot sampling hint from the p.d.f. that defines the observable // obs - for(auto const& func : _actualVars) { - auto hint = static_cast(func)->plotSamplingHint(obs, xlo, xhi); + for(auto *func : static_range_cast(_actualVars)) { + auto hint = func->plotSamplingHint(obs, xlo, xhi); if (hint) { return hint; } diff --git a/roofit/roofitcore/src/RooSuperCategory.cxx b/roofit/roofitcore/src/RooSuperCategory.cxx index 00f20eb86d6aa..e8fcbe5df667a 100644 --- a/roofit/roofitcore/src/RooSuperCategory.cxx +++ b/roofit/roofitcore/src/RooSuperCategory.cxx @@ -93,8 +93,7 @@ bool RooSuperCategory::setIndex(Int_t index, bool printError) } bool error = false; - for (auto arg : _multiCat->_catSet) { - auto cat = static_cast(arg); + for (auto* cat : static_range_cast(_multiCat->_catSet)) { if (cat->empty()) { if (printError) { coutE(InputArguments) << __func__ << ": Found a category with zero states. Cannot set state for '" diff --git a/roofit/roostats/src/BayesianCalculator.cxx b/roofit/roostats/src/BayesianCalculator.cxx index 3d5d2e513231d..c6602e432e8bb 100644 --- a/roofit/roostats/src/BayesianCalculator.cxx +++ b/roofit/roostats/src/BayesianCalculator.cxx @@ -833,8 +833,7 @@ RooAbsReal* BayesianCalculator::GetPosteriorFunction() const << " Negative log likelihood evaluates to infinity " << std::endl << " Non-const Parameter values : "; RooArgList p(*constrainedParams); - for (std::size_t i = 0; i < p.size(); ++i) { - RooRealVar * v = dynamic_cast(&p[i] ); + for (auto *v : dynamic_range_cast(p)) { if (v!=nullptr) ccoutE(Eval) << v->GetName() << " = " << v->getVal() << " "; } ccoutE(Eval) << std::endl; diff --git a/roofit/roostats/src/DetailedOutputAggregator.cxx b/roofit/roostats/src/DetailedOutputAggregator.cxx index 6686f536b689b..8868038fa1571 100644 --- a/roofit/roostats/src/DetailedOutputAggregator.cxx +++ b/roofit/roostats/src/DetailedOutputAggregator.cxx @@ -131,8 +131,8 @@ namespace RooStats { } fResult->add(RooArgSet(*fBuiltSet), weight); - for (RooAbsArg* v : *fBuiltSet) { - if (RooRealVar* var= dynamic_cast(v)) { + for (auto* var : dynamic_range_cast(*fBuiltSet)) { + if (var) { // Invalidate values in case we don't set some of them next time round (eg. if fit not done) var->setVal(std::numeric_limits::quiet_NaN()); var->removeError(); diff --git a/roofit/roostats/src/HypoTestInverter.cxx b/roofit/roostats/src/HypoTestInverter.cxx index 16f12bb86d21a..683f913b81750 100644 --- a/roofit/roostats/src/HypoTestInverter.cxx +++ b/roofit/roostats/src/HypoTestInverter.cxx @@ -1225,8 +1225,7 @@ SamplingDistribution * HypoTestInverter::RebuildDistributions(bool isUpper, int RooArgList genObs(*bkgdata->get(0)); RooStats::PrintListContent(genObs, oocoutP(nullptr,Generation) ); nObs = 0; - for (std::size_t i = 0; i < genObs.size(); ++i) { - RooRealVar * x = dynamic_cast(&genObs[i]); + for (auto *x : dynamic_range_cast(genObs)) { if (x) nObs += x->getVal(); } } diff --git a/roofit/roostats/src/LikelihoodInterval.cxx b/roofit/roostats/src/LikelihoodInterval.cxx index 01089ae208fd5..0ec1e5e3fb19e 100644 --- a/roofit/roostats/src/LikelihoodInterval.cxx +++ b/roofit/roostats/src/LikelihoodInterval.cxx @@ -233,12 +233,11 @@ bool LikelihoodInterval::CreateMinimizer() { // need to restore values and errors for POI if (fBestFitParams) { - for (std::size_t i = 0; i < params.size(); ++i) { - RooRealVar & par = static_cast( params[i]); - RooRealVar * fitPar = static_cast (fBestFitParams->find(par.GetName() ) ); + for (auto *par : static_range_cast(params)) { + RooRealVar * fitPar = static_cast (fBestFitParams->find(par->GetName() ) ); if (fitPar) { - par.setVal( fitPar->getVal() ); - par.setError( fitPar->getError() ); + par->setVal( fitPar->getVal() ); + par->setError( fitPar->getError() ); } } } diff --git a/roofit/roostats/src/LikelihoodIntervalPlot.cxx b/roofit/roostats/src/LikelihoodIntervalPlot.cxx index d072e35d1459a..d21e5e2b867b4 100644 --- a/roofit/roostats/src/LikelihoodIntervalPlot.cxx +++ b/roofit/roostats/src/LikelihoodIntervalPlot.cxx @@ -340,11 +340,10 @@ void LikelihoodIntervalPlot::Draw(const Option_t *options) std::unique_ptr vars{newProfile->getVariables()}; RooArgList params(*vars); // set values and error for the POI to the best fit values - for (std::size_t i = 0; i < params.size(); ++i) { - RooRealVar & par = static_cast( params[i]); - RooRealVar * fitPar = static_cast (fInterval->GetBestFitParameters()->find(par.GetName() ) ); + for (auto *par : static_range_cast(params)) { + RooRealVar * fitPar = static_cast (fInterval->GetBestFitParameters()->find(par->GetName() ) ); if (fitPar) { - par.setVal( fitPar->getVal() ); + par->setVal( fitPar->getVal() ); } } // do a profile evaluation to start from the best fit values of parameters diff --git a/roofit/roostats/src/ProfileLikelihoodCalculator.cxx b/roofit/roostats/src/ProfileLikelihoodCalculator.cxx index 0faadf305b13f..8ceced2cb0c41 100644 --- a/roofit/roostats/src/ProfileLikelihoodCalculator.cxx +++ b/roofit/roostats/src/ProfileLikelihoodCalculator.cxx @@ -237,12 +237,11 @@ LikelihoodInterval* ProfileLikelihoodCalculator::GetInterval() const { // t.b.f. " RooProfileLL should keep and provide possibility to query on global minimum // set POI to fit value (this will speed up profileLL calculation of global minimum) const RooArgList & fitParams = fFitResult->floatParsFinal(); - for (std::size_t i = 0; i < fitParams.size(); ++i) { - RooRealVar & fitPar = static_cast( fitParams[i]); - RooRealVar * par = static_cast(fPOI.find( fitPar.GetName() )); + for (auto *fitPar : static_range_cast(fitParams)) { + RooRealVar * par = static_cast(fPOI.find( fitPar->GetName() )); if (par) { - par->setVal( fitPar.getVal() ); - par->setError( fitPar.getError() ); + par->setVal( fitPar->getVal() ); + par->setError( fitPar->getError() ); } } From 69b34b7cc09f9a913541926129d65e87a1e9288b Mon Sep 17 00:00:00 2001 From: Carsten Burgard Date: Wed, 17 Jun 2026 17:36:07 +0200 Subject: [PATCH 4/5] [RF][HS3] Added RooWrapperPdf and updated extendpdf to be HS3-compliant Added I/O for RooWrapperPdf and changed extendpdf I/O to be compliant with HS3. (cherry picked from commit d75a72f9b3ef9db478f7159e8f01cefada212e88) --- roofit/hs3/src/JSONFactories_RooFitCore.cxx | 54 ++++++++++++++++++++- 1 file changed, 53 insertions(+), 1 deletion(-) diff --git a/roofit/hs3/src/JSONFactories_RooFitCore.cxx b/roofit/hs3/src/JSONFactories_RooFitCore.cxx index 29a3f66fa99d4..7fae03dc51b2b 100644 --- a/roofit/hs3/src/JSONFactories_RooFitCore.cxx +++ b/roofit/hs3/src/JSONFactories_RooFitCore.cxx @@ -48,6 +48,7 @@ #include #include #include +#include #include #include #include @@ -1187,17 +1188,65 @@ class RooSplineStreamer : public RooFit::JSONIO::Exporter { } }; +class RooWrapperPdfImporter : public RooFit::JSONIO::Importer { +public: + bool importArg(RooJSONFactoryWSTool *tool, const RooFit::Detail::JSONNode &node) const override + { + if (node["type"].val() != "density_function_dist") + return false; + + auto name = RooJSONFactoryWSTool::name(node); + auto *func = tool->requestArg(node, "function"); + + bool selfNormalized = false; + + if (auto sn = node.find("selfNormalized")) + selfNormalized = sn->val_bool(); + + tool->wsEmplace(name, *func, selfNormalized); + return true; + } +}; + +class RooWrapperPdfStreamer : public RooFit::JSONIO::Exporter { +public: + std::string const &key() const override; + + bool exportObject(RooJSONFactoryWSTool *, const RooAbsArg *arg, RooFit::Detail::JSONNode &node) const override + { + auto const *pdf = dynamic_cast(arg); + if (!pdf) + return false; + + node["type"] << "density_function_dist"; + + // Proxy name in RooWrapperPdf is "_func" / "func" depending on accessor/proxy export. + // Prefer a public accessor if one exists; otherwise inspect proxies as below. + auto const *funcProxy = dynamic_cast(pdf->getProxy(0)); + if (!funcProxy || !funcProxy->absArg()) + return false; + + node["function"] << funcProxy->absArg()->GetName(); + if (pdf->selfNormalized()) + node["selfnormalized"] << true; + + return true; + } +}; + #define DEFINE_EXPORTER_KEY(class_name, name) \ std::string const &class_name::key() const \ { \ const static std::string keystring = name; \ return keystring; \ } + template <> DEFINE_EXPORTER_KEY(RooAddPdfStreamer, "mixture_dist"); template <> DEFINE_EXPORTER_KEY(RooAddPdfStreamer, "mixture_model"); DEFINE_EXPORTER_KEY(RooBinSamplingPdfStreamer, "binsampling"); +DEFINE_EXPORTER_KEY(RooWrapperPdfStreamer, "density_function_dist"); DEFINE_EXPORTER_KEY(RooBinWidthFunctionStreamer, "binwidth"); DEFINE_EXPORTER_KEY(RooLegacyExpPolyStreamer, "legacy_exp_poly_dist"); DEFINE_EXPORTER_KEY(RooExponentialStreamer, "exponential_dist"); @@ -1223,7 +1272,7 @@ DEFINE_EXPORTER_KEY(RooTFnBindingStreamer, "generic_function"); DEFINE_EXPORTER_KEY(RooRealIntegralStreamer, "integral"); DEFINE_EXPORTER_KEY(RooDerivativeStreamer, "derivative"); DEFINE_EXPORTER_KEY(RooFFTConvPdfStreamer, "fft_conv_pdf"); -DEFINE_EXPORTER_KEY(RooExtendPdfStreamer, "extend_pdf"); +DEFINE_EXPORTER_KEY(RooExtendPdfStreamer, "rate_extended_dist"); DEFINE_EXPORTER_KEY(ParamHistFuncStreamer, "step"); DEFINE_EXPORTER_KEY(RooSplineStreamer, "spline"); @@ -1234,6 +1283,8 @@ DEFINE_EXPORTER_KEY(RooSplineStreamer, "spline"); STATIC_EXECUTE([]() { using namespace RooFit::JSONIO; + registerImporter("density_function_dist"); + registerImporter("rate_extended_dist"); registerImporter("product", false); registerImporter("product_dist", false); registerImporter("sum", false); @@ -1264,6 +1315,7 @@ STATIC_EXECUTE([]() { registerImporter("step", false); registerImporter("spline", false); + registerExporter(RooWrapperPdf::Class()); registerExporter>(RooAddPdf::Class(), false); registerExporter>(RooAddModel::Class(), false); registerExporter(RooBinSamplingPdf::Class(), false); From d8bbec359712aa839278135bceef0b6cb6b25c7e Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Fri, 19 Jun 2026 20:28:18 +0200 Subject: [PATCH 5/5] [RF][HS3] Reduce code duplication in RooFit HS3 Refactor the HS3 serialization code to remove duplicated logic without any change in behavior. All hs3, histfactory and JSON tests still pass. JSONFactories_RooFitCore.cxx: - Drop RooLegacyExpPolyFactory; it was identical to the existing RooPolynomialFactory, which is now registered directly. - Drop RooLegacyExpPolyStreamer; it was identical to RooPolynomialStreamer. - Merge RooHistFuncStreamer/RooHistPdfStreamer into RooHistStreamer and RooHistFuncFactory/RooHistPdfFactory into RooHistFactory; each pair differed only by the concrete RooFit type. RooJSONFactoryWSTool.cxx: - Remove fillSeqSanitizedName, which had no callers anywhere in ROOT. - Factor the "make a valid name but refuse to change the first character" block (exportCategory, exportCombinedData) into makeValidNameOrError(). - Collapse the repeated ofstream/ifstream open-and-error boilerplate in exportJSON/exportYML/importJSON/importYML. error() is [[noreturn]], so the trailing 'return false' statements were dead code. JSONFactories_HistFactory.cxx: - Remove a redundant mid-file block of #includes ( was already included; the others were already used earlier in the translation unit). JSONIO.cxx: - Deduplicate removeImporters/removeExporters and printImporters/ printExporters into shared templated helpers. (cherry picked from commit b763ff1d1bb3325f3d981ff96c90a48517f6c6fa) --- .../hs3/inc/RooFitHS3/RooJSONFactoryWSTool.h | 1 - roofit/hs3/src/JSONFactories_HistFactory.cxx | 6 -- roofit/hs3/src/JSONFactories_RooFitCore.cxx | 101 ++++-------------- roofit/hs3/src/JSONIO.cxx | 53 +++++---- roofit/hs3/src/RooJSONFactoryWSTool.cxx | 97 +++++------------ 5 files changed, 72 insertions(+), 186 deletions(-) diff --git a/roofit/hs3/inc/RooFitHS3/RooJSONFactoryWSTool.h b/roofit/hs3/inc/RooFitHS3/RooJSONFactoryWSTool.h index c58ce70154cfc..9fdb0c530d53c 100644 --- a/roofit/hs3/inc/RooFitHS3/RooJSONFactoryWSTool.h +++ b/roofit/hs3/inc/RooFitHS3/RooJSONFactoryWSTool.h @@ -63,7 +63,6 @@ class RooJSONFactoryWSTool { static RooFit::Detail::JSONNode const *findNamedChild(RooFit::Detail::JSONNode const &node, std::string const &name); static void fillSeq(RooFit::Detail::JSONNode &node, RooAbsCollection const &coll, size_t nMax = -1); - static void fillSeqSanitizedName(RooFit::Detail::JSONNode &node, RooAbsCollection const &coll, size_t nMax = -1); template T *request(const std::string &objname, const std::string &requestAuthor) diff --git a/roofit/hs3/src/JSONFactories_HistFactory.cxx b/roofit/hs3/src/JSONFactories_HistFactory.cxx index f3cffdbb7a9f5..ea1b3de1672af 100644 --- a/roofit/hs3/src/JSONFactories_HistFactory.cxx +++ b/roofit/hs3/src/JSONFactories_HistFactory.cxx @@ -765,12 +765,6 @@ std::vector splitTopLevelProduct(const std::string &expr) return parts; } -#include -#include -#include -#include -#include - NormSys parseOverallModifierFormula(const std::string &s, RooFormulaVar *formula) { static const std::regex pattern( diff --git a/roofit/hs3/src/JSONFactories_RooFitCore.cxx b/roofit/hs3/src/JSONFactories_RooFitCore.cxx index 7fae03dc51b2b..42fe920e04e18 100644 --- a/roofit/hs3/src/JSONFactories_RooFitCore.cxx +++ b/roofit/hs3/src/JSONFactories_RooFitCore.cxx @@ -499,37 +499,6 @@ class RooExponentialFactory : public RooFit::JSONIO::Importer { } }; -class RooLegacyExpPolyFactory : public RooFit::JSONIO::Importer { -public: - bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override - { - std::string name(RooJSONFactoryWSTool::name(p)); - if (!p.has_child("coefficients")) { - RooJSONFactoryWSTool::error("no coefficients given in '" + name + "'"); - } - RooAbsReal *x = tool->requestArg(p, "x"); - RooArgList coefs; - int order = 0; - int lowestOrder = 0; - for (const auto &coef : p["coefficients"].children()) { - // As long as the coefficients match the default coefficients in - // RooFit, we don't have to instantiate RooFit objects but can - // increase the lowestOrder flag. - if (order == 0 && coef.val() == "1.0") { - ++lowestOrder; - } else if (coefs.empty() && coef.val() == "0.0") { - ++lowestOrder; - } else { - coefs.add(*tool->request(coef.val(), name)); - } - ++order; - } - - tool->wsEmplace(name, *x, coefs, lowestOrder); - return true; - } -}; - class RooMultiVarGaussianFactory : public RooFit::JSONIO::Importer { public: bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override @@ -743,40 +712,13 @@ class RooRealSumFuncStreamer : public RooFit::JSONIO::Exporter { } }; -class RooHistFuncStreamer : public RooFit::JSONIO::Exporter { -public: - std::string const &key() const override; - bool exportObject(RooJSONFactoryWSTool *tool, const RooAbsArg *func, JSONNode &elem) const override - { - const RooHistFunc *hf = static_cast(func); - elem["type"] << key(); - RooDataHist const &dh = hf->dataHist(); - tool->exportHisto(*dh.get(), dh.numEntries(), dh.weightArray(), elem["data"].set_map()); - return true; - } -}; - -class RooHistFuncFactory : public RooFit::JSONIO::Importer { -public: - bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override - { - std::string name(RooJSONFactoryWSTool::name(p)); - if (!p.has_child("data")) { - RooJSONFactoryWSTool::error("function '" + name + "' is of histogram type, but does not define a 'data' key"); - } - std::unique_ptr dataHist = - RooJSONFactoryWSTool::readBinnedData(p["data"], name, RooJSONFactoryWSTool::readAxes(p["data"])); - tool->wsEmplace(name, *dataHist->get(), *dataHist); - return true; - } -}; - -class RooHistPdfStreamer : public RooFit::JSONIO::Exporter { +template +class RooHistStreamer : public RooFit::JSONIO::Exporter { public: std::string const &key() const override; bool exportObject(RooJSONFactoryWSTool *tool, const RooAbsArg *func, JSONNode &elem) const override { - const RooHistPdf *hf = static_cast(func); + const RooArg_t *hf = static_cast(func); elem["type"] << key(); RooDataHist const &dh = hf->dataHist(); tool->exportHisto(*dh.get(), dh.numEntries(), dh.weightArray(), elem["data"].set_map()); @@ -784,7 +726,8 @@ class RooHistPdfStreamer : public RooFit::JSONIO::Exporter { } }; -class RooHistPdfFactory : public RooFit::JSONIO::Importer { +template +class RooHistFactory : public RooFit::JSONIO::Importer { public: bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override { @@ -794,7 +737,7 @@ class RooHistPdfFactory : public RooFit::JSONIO::Importer { } std::unique_ptr dataHist = RooJSONFactoryWSTool::readBinnedData(p["data"], name, RooJSONFactoryWSTool::readAxes(p["data"])); - tool->wsEmplace(name, *dataHist->get(), *dataHist); + tool->wsEmplace(name, *dataHist->get(), *dataHist); return true; } }; @@ -893,17 +836,6 @@ class RooPolynomialStreamer : public RooFit::JSONIO::Exporter { } }; -class RooLegacyExpPolyStreamer : public RooFit::JSONIO::Exporter { -public: - std::string const &key() const override; - bool exportObject(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &elem) const override - { - elem["type"] << key(); - writePolynomialBody(static_cast(func), elem); - return true; - } -}; - class RooPoissonStreamer : public RooFit::JSONIO::Exporter { public: std::string const &key() const override; @@ -1248,14 +1180,17 @@ DEFINE_EXPORTER_KEY(RooAddPdfStreamer, "mixture_model"); DEFINE_EXPORTER_KEY(RooBinSamplingPdfStreamer, "binsampling"); DEFINE_EXPORTER_KEY(RooWrapperPdfStreamer, "density_function_dist"); DEFINE_EXPORTER_KEY(RooBinWidthFunctionStreamer, "binwidth"); -DEFINE_EXPORTER_KEY(RooLegacyExpPolyStreamer, "legacy_exp_poly_dist"); +template <> +DEFINE_EXPORTER_KEY(RooPolynomialStreamer, "legacy_exp_poly_dist"); DEFINE_EXPORTER_KEY(RooExponentialStreamer, "exponential_dist"); template <> DEFINE_EXPORTER_KEY(RooFormulaArgStreamer, "generic_function"); template <> DEFINE_EXPORTER_KEY(RooFormulaArgStreamer, "generic_dist"); -DEFINE_EXPORTER_KEY(RooHistFuncStreamer, "histogram"); -DEFINE_EXPORTER_KEY(RooHistPdfStreamer, "histogram_dist"); +template <> +DEFINE_EXPORTER_KEY(RooHistStreamer, "histogram"); +template <> +DEFINE_EXPORTER_KEY(RooHistStreamer, "histogram_dist"); DEFINE_EXPORTER_KEY(RooLogNormalStreamer, "lognormal_dist"); DEFINE_EXPORTER_KEY(RooMultiVarGaussianStreamer, "multivariate_normal_dist"); DEFINE_EXPORTER_KEY(RooPoissonStreamer, "poisson_dist"); @@ -1292,12 +1227,12 @@ STATIC_EXECUTE([]() { registerImporter("mixture_model", false); registerImporter("binsampling_dist", false); registerImporter("binwidth", false); - registerImporter("legacy_exp_poly_dist", false); + registerImporter>("legacy_exp_poly_dist", false); registerImporter("exponential_dist", false); registerImporter>("generic_function", false); registerImporter>("generic_dist", false); - registerImporter("histogram", false); - registerImporter("histogram_dist", false); + registerImporter>("histogram", false); + registerImporter>("histogram_dist", false); registerImporter("lognormal_dist", false); registerImporter("multivariate_normal_dist", false); registerImporter("poisson_dist", false); @@ -1320,12 +1255,12 @@ STATIC_EXECUTE([]() { registerExporter>(RooAddModel::Class(), false); registerExporter(RooBinSamplingPdf::Class(), false); registerExporter(RooBinWidthFunction::Class(), false); - registerExporter(RooLegacyExpPoly::Class(), false); + registerExporter>(RooLegacyExpPoly::Class(), false); registerExporter(RooExponential::Class(), false); registerExporter>(RooFormulaVar::Class(), false); registerExporter>(RooGenericPdf::Class(), false); - registerExporter(RooHistFunc::Class(), false); - registerExporter(RooHistPdf::Class(), false); + registerExporter>(RooHistFunc::Class(), false); + registerExporter>(RooHistPdf::Class(), false); registerExporter(RooLognormal::Class(), false); registerExporter(RooMultiVarGaussian::Class(), false); registerExporter(RooPoisson::Class(), false); diff --git a/roofit/hs3/src/JSONIO.cxx b/roofit/hs3/src/JSONIO.cxx index 5381887730a17..8d6cc4099624d 100644 --- a/roofit/hs3/src/JSONIO.cxx +++ b/roofit/hs3/src/JSONIO.cxx @@ -147,26 +147,13 @@ bool registerExporter(const std::string &key, std::unique_ptr f, return true; } -int removeImporters(const std::string &needle) -{ - int n = 0; - for (auto &element : importers()) { - for (size_t i = element.second.size(); i > 0; --i) { - auto *imp = element.second[i - 1].get(); - std::string name(typeid(*imp).name()); - if (name.find(needle) != std::string::npos) { - element.second.erase(element.second.begin() + i - 1); - ++n; - } - } - } - return n; -} +namespace { -int removeExporters(const std::string &needle) +template +int removeByTypeName(Map &map, const std::string &needle) { int n = 0; - for (auto &element : exporters()) { + for (auto &element : map) { for (size_t i = element.second.size(); i > 0; --i) { auto *imp = element.second[i - 1].get(); std::string name(typeid(*imp).name()); @@ -179,25 +166,37 @@ int removeExporters(const std::string &needle) return n; } -void printImporters() +template +void printByTypeName(Map &map, KeyPrinter printKey) { - for (const auto &x : importers()) { + for (const auto &x : map) { for (const auto &ePtr : x.second) { // Passing *e directory to typeid results in clang warnings. auto const &e = *ePtr; - std::cout << x.first << "\t" << typeid(e).name() << std::endl; + std::cout << printKey(x.first) << "\t" << typeid(e).name() << std::endl; } } } + +} // namespace + +int removeImporters(const std::string &needle) +{ + return removeByTypeName(importers(), needle); +} + +int removeExporters(const std::string &needle) +{ + return removeByTypeName(exporters(), needle); +} + +void printImporters() +{ + printByTypeName(importers(), [](auto const &key) { return key; }); +} void printExporters() { - for (const auto &x : exporters()) { - for (const auto &ePtr : x.second) { - // Passing *e directory to typeid results in clang warnings. - auto const &e = *ePtr; - std::cout << x.first->GetName() << "\t" << typeid(e).name() << std::endl; - } - } + printByTypeName(exporters(), [](auto const &key) { return key->GetName(); }); } void loadFactoryExpressions(const std::string &fname) diff --git a/roofit/hs3/src/RooJSONFactoryWSTool.cxx b/roofit/hs3/src/RooJSONFactoryWSTool.cxx index 9d7acb2e0c94d..65634d6da2a40 100644 --- a/roofit/hs3/src/RooJSONFactoryWSTool.cxx +++ b/roofit/hs3/src/RooJSONFactoryWSTool.cxx @@ -909,26 +909,6 @@ void RooJSONFactoryWSTool::fillSeq(JSONNode &node, RooAbsCollection const &coll, } } -void RooJSONFactoryWSTool::fillSeqSanitizedName(JSONNode &node, RooAbsCollection const &coll, size_t nMax) -{ - const size_t old_children = node.num_children(); - node.set_seq(); - size_t n = 0; - for (RooAbsArg const *arg : coll) { - if (n >= nMax) - break; - if (isLiteralConstVar(*arg)) { - node.append_child() << static_cast(arg)->getVal(); - } else { - node.append_child() << sanitizeName(arg->GetName()); - } - ++n; - } - if (node.num_children() != old_children + coll.size()) { - error("unable to stream collection " + std::string(coll.GetName()) + " to " + node.key()); - } -} - JSONNode &RooJSONFactoryWSTool::appendNamedChild(JSONNode &node, std::string const &name) { if (!useListsInsteadOfDicts) { @@ -1489,25 +1469,31 @@ void RooJSONFactoryWSTool::exportArray(std::size_t n, double const *contents, JS * @param node The JSONNode to which the category data will be exported. * @return void */ +namespace { + +// Turn an arbitrary string into a valid variable name, but refuse to change the +// first character (which would silently rename the object). +std::string makeValidNameOrError(std::string const &in) +{ + if (!std::isalpha(in[0])) { + RooJSONFactoryWSTool::error("refusing to change first character of string '" + in + "' to make a valid name!"); + } + std::string out = RooFit::Detail::makeValidVarName(in); + if (out != in) { + oocoutW(nullptr, IO) << "RooFitHS3: changed '" << in << "' to '" << out << "' to become a valid name"; + } + return out; +} + +} // namespace + void RooJSONFactoryWSTool::exportCategory(RooAbsCategory const &cat, JSONNode &node) { auto &labels = node["labels"].set_seq(); auto &indices = node["indices"].set_seq(); for (auto const &item : cat) { - std::string label; - if (std::isalpha(item.first[0])) { - label = RooFit::Detail::makeValidVarName(item.first); - if (label != item.first) { - oocoutW(nullptr, IO) << "RooFitHS3: changed '" << item.first << "' to '" << label - << "' to become a valid name"; - } - } else { - RooJSONFactoryWSTool::error("refusing to change first character of string '" + item.first + - "' to make a valid name!"); - label = item.first; - } - labels.append_child() << label; + labels.append_child() << makeValidNameOrError(item.first); indices.append_child() << item.second; } } @@ -1569,18 +1555,7 @@ RooJSONFactoryWSTool::CombinedData RooJSONFactoryWSTool::exportCombinedData(RooA for (std::unique_ptr const &absData : dataList) { std::string catName(absData->GetName()); - std::string dataName; - if (std::isalpha(catName[0])) { - dataName = RooFit::Detail::makeValidVarName(catName); - if (dataName != catName) { - oocoutW(nullptr, IO) << "RooFitHS3: changed '" << catName << "' to '" << dataName - << "' to become a valid name"; - } - } else { - RooJSONFactoryWSTool::error("refusing to change first character of string '" + catName + - "' to make a valid name!"); - dataName = catName; - } + std::string dataName = makeValidNameOrError(catName); absData->SetName((std::string(data.GetName()) + "_" + dataName).c_str()); datamap.components[catName] = absData->GetName(); this->exportData(*absData); @@ -2162,12 +2137,8 @@ bool RooJSONFactoryWSTool::exportJSON(std::ostream &os) bool RooJSONFactoryWSTool::exportJSON(std::string const &filename) { std::ofstream out(filename.c_str()); - if (!out.is_open()) { - std::stringstream ss; - ss << "RooJSONFactoryWSTool() invalid output file '" << filename << "'." << std::endl; - RooJSONFactoryWSTool::error(ss.str()); - return false; - } + if (!out.is_open()) + RooJSONFactoryWSTool::error("RooJSONFactoryWSTool() invalid output file '" + filename + "'."); return this->exportJSON(out); } @@ -2195,12 +2166,8 @@ bool RooJSONFactoryWSTool::exportYML(std::ostream &os) bool RooJSONFactoryWSTool::exportYML(std::string const &filename) { std::ofstream out(filename.c_str()); - if (!out.is_open()) { - std::stringstream ss; - ss << "RooJSONFactoryWSTool() invalid output file '" << filename << "'." << std::endl; - RooJSONFactoryWSTool::error(ss.str()); - return false; - } + if (!out.is_open()) + RooJSONFactoryWSTool::error("RooJSONFactoryWSTool() invalid output file '" + filename + "'."); return this->exportYML(out); } @@ -2405,12 +2372,8 @@ bool RooJSONFactoryWSTool::importJSON(std::string const &filename) { // import a JSON file to the workspace std::ifstream infile(filename.c_str()); - if (!infile.is_open()) { - std::stringstream ss; - ss << "RooJSONFactoryWSTool() invalid input file '" << filename << "'." << std::endl; - RooJSONFactoryWSTool::error(ss.str()); - return false; - } + if (!infile.is_open()) + RooJSONFactoryWSTool::error("RooJSONFactoryWSTool() invalid input file '" + filename + "'."); return this->importJSON(infile); } @@ -2440,12 +2403,8 @@ bool RooJSONFactoryWSTool::importYML(std::string const &filename) { // import a YML file to the workspace std::ifstream infile(filename.c_str()); - if (!infile.is_open()) { - std::stringstream ss; - ss << "RooJSONFactoryWSTool() invalid input file '" << filename << "'." << std::endl; - RooJSONFactoryWSTool::error(ss.str()); - return false; - } + if (!infile.is_open()) + RooJSONFactoryWSTool::error("RooJSONFactoryWSTool() invalid input file '" + filename + "'."); return this->importYML(infile); }