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/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 df59000016fc0..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( @@ -1412,8 +1406,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..42fe920e04e18 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 @@ -498,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 @@ -627,8 +597,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; @@ -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; @@ -1188,26 +1120,77 @@ 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"); +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"); @@ -1224,7 +1207,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"); @@ -1235,6 +1218,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); @@ -1242,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); @@ -1265,16 +1250,17 @@ STATIC_EXECUTE([]() { registerImporter("step", false); registerImporter("spline", false); + registerExporter(RooWrapperPdf::Class()); registerExporter>(RooAddPdf::Class(), false); 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); } 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/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/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/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/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/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/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/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/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 9fa2bed777aaf..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(); } } @@ -1237,7 +1236,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/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 5d4d669f5aa28..d21e5e2b867b4 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,13 +337,13 @@ 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]); - 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/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/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() ); } } 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/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 @@ -8905,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"); } @@ -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