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