diff --git a/roofit/histfactory/test/testHistFactory.cxx b/roofit/histfactory/test/testHistFactory.cxx index 29acd2c4ded64..252ec9e7fed50 100644 --- a/roofit/histfactory/test/testHistFactory.cxx +++ b/roofit/histfactory/test/testHistFactory.cxx @@ -9,6 +9,7 @@ #endif #include +#include #include #include @@ -18,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -32,7 +34,11 @@ #include "../../roofitcore/test/gtest_wrapper.h" +#include +#include +#include #include +#include #include namespace { @@ -138,6 +144,105 @@ TEST(HistFactory, Read_ROOT6_16_Combined_Model) EXPECT_NEAR(pdf->getVal(*obs), 0.95652174, 1.E-8); } +// The reference workspace round-trips through HS3 with its constraints collected into array distributions: the +// re-export gives the same document, and the likelihood is the same. +TEST(HistFactory, HS3ArrayDistRoundTrip_ROOT6_16_Combined_Model) +{ + RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING, 0u, RooFit::NumericIntegration, true}; + + std::string filename = "./ref_6.16_example_UsingC_combined_meas_model.root"; + std::unique_ptr file(TFile::Open(filename.c_str())); + ASSERT_TRUE(file && file->IsOpen()); + RooWorkspace *ws; + file->GetObject("combined", ws); + ASSERT_NE(ws, nullptr); + + const std::string js = RooJSONFactoryWSTool{*ws}.exportJSONtoString(); + auto tree = RooFit::Detail::JSONTree::create(js); + auto const &distributions = tree->rootnode()["distributions"]; + for (auto const &dist : distributions.children()) { + EXPECT_NE(dist["type"].val(), "gaussian_dist") << js; + EXPECT_NE(dist["type"].val(), "poisson_dist") << js; + } + const auto *gauss = RooJSONFactoryWSTool::findNamedChild(distributions, "simPdf_gaussian_constraints"); + ASSERT_NE(gauss, nullptr) << js; + EXPECT_EQ((*gauss)["components"].num_children(), 4u) << js; // Lumi, alpha_syst1, alpha_syst2, alpha_syst3 + const auto *pois = RooJSONFactoryWSTool::findNamedChild(distributions, "simPdf_poisson_constraints"); + ASSERT_NE(pois, nullptr) << js; + EXPECT_EQ((*pois)["components"].num_children(), 2u) << js; + // HistFactory builds the Gaussian constraints with the parameter as x and the global observable as mean; they + // are written the other way round, so that the observable is x. Also for the constant parameter Lumi. + std::map gaussMeans; + for (auto const &component : (*gauss)["components"].children()) { + gaussMeans[component["x"].val()] = component["mean"].val(); + } + EXPECT_EQ(gaussMeans, (std::map{{"nominalLumi", "Lumi"}, + {"nom_alpha_syst1", "alpha_syst1"}, + {"nom_alpha_syst2", "alpha_syst2"}, + {"nom_alpha_syst3", "alpha_syst3"}})) + << js; + + // The parameter sections of this workspace from ROOT 6.16 are not stable under a round trip independently of the + // arrays (its legacy bin-width variables come back as different objects), so the model sections are compared. + RooWorkspace wsFromJson("new"); + RooJSONFactoryWSTool{wsFromJson}.importJSONfromString(js); + const std::string js2 = RooJSONFactoryWSTool{wsFromJson}.exportJSONtoString(); + auto section = [](std::string const &json, const char *key) { + auto t = RooFit::Detail::JSONTree::create(json); + std::stringstream ss; + t->rootnode()[key].writeJSON(ss); + return ss.str(); + }; + for (const char *key : {"distributions", "functions", "likelihoods", "analyses"}) { + EXPECT_EQ(section(js2, key), section(js, key)) << key; + } + // The re-imported Gaussians have the observable as x. + auto *reimported = dynamic_cast(wsFromJson.pdf("nom_alpha_syst2_constraint")); + ASSERT_NE(reimported, nullptr); + EXPECT_STREQ(reimported->getX().GetName(), "nom_alpha_syst2"); + EXPECT_STREQ(reimported->getMean().GetName(), "alpha_syst2"); + auto *lumi = dynamic_cast(wsFromJson.pdf("nominalLumi_constraint")); + ASSERT_NE(lumi, nullptr); + EXPECT_STREQ(lumi->getX().GetName(), "nominalLumi"); + EXPECT_STREQ(lumi->getMean().GetName(), "Lumi"); + + // The same without arrays: the standalone gaussian_dist nodes have the observable as x too. + RooJSONFactoryWSTool::config().exportArrayDist = false; + const std::string plain = RooJSONFactoryWSTool{*ws}.exportJSONtoString(); + RooJSONFactoryWSTool::config().exportArrayDist = true; + auto plainTree = RooFit::Detail::JSONTree::create(plain); + const auto *standalone = + RooJSONFactoryWSTool::findNamedChild(plainTree->rootnode()["distributions"], "alpha_syst2Constraint"); + ASSERT_NE(standalone, nullptr) << plain; + EXPECT_EQ((*standalone)["type"].val(), "gaussian_dist") << plain; + EXPECT_EQ((*standalone)["x"].val(), "nom_alpha_syst2") << plain; + EXPECT_EQ((*standalone)["mean"].val(), "alpha_syst2") << plain; + RooWorkspace wsPlain("plain"); + RooJSONFactoryWSTool{wsPlain}.importJSONfromString(plain); + auto *plainGauss = dynamic_cast(wsPlain.pdf("alpha_syst2Constraint")); + ASSERT_NE(plainGauss, nullptr); + EXPECT_STREQ(plainGauss->getX().GetName(), "nom_alpha_syst2"); + + // The likelihood is compared up to a constant: whether RooFit uses its binned likelihood for the channels (which + // differs from the event-wise one by a constant) is not part of the document. + auto nllDifference = [](RooWorkspace &w) { + auto *mc = dynamic_cast(w.obj("ModelConfig")); + if (!mc || !w.data("obsData")) { + return std::numeric_limits::quiet_NaN(); + } + std::unique_ptr nll{ + mc->GetPdf()->createNLL(*w.data("obsData"), RooFit::GlobalObservables(*mc->GetGlobalObservables()))}; + const double nominal = nll->getVal(); + w.var("alpha_syst2")->setVal(0.5); + w.var("gamma_stat_channel1_bin_0")->setVal(1.1); + return nll->getVal() - nominal; + }; + const double reference = nllDifference(*ws); + EXPECT_FALSE(std::isnan(reference)); + EXPECT_NE(reference, 0.0); + EXPECT_NEAR(nllDifference(wsFromJson), reference, 1e-8 * std::max(1.0, std::abs(reference))); +} + /// What kind of model is set up. Use this to instantiate /// a test suite. enum class MakeModelMode { @@ -700,15 +805,10 @@ INSTANTIATE_TEST_SUITE_P(HistFactory, HFFixtureFit, testing::Values(ROOFIT_EVAL_BACKENDS_WITH_CODEGEN)), getNameFromInfo); -// Regression test for the HS3 importer's handling of "shapefactor" modifiers. -// HistFactory's ShapeFactor is exported with type "shapesys" (the only -// modifier type the writer ever emits for a ParamHistFunc), but valid HS3 JSON -// may also use the dedicated "shapefactor" type for an unconstrained -// ParamHistFunc - and the importer must accept it. This test starts from a -// HistFactory model built with MakeModelAndMeasurementFast, rewrites the -// modifier type to "shapefactor" in the JSON, and checks that the importer -// recognises it. Without the fix in JSONFactories_HistFactory.cxx the import -// throws "modifier ... of unknown type 'shapefactor'". +// HS3 round trip of a free HistFactory ShapeFactor. It is written as a +// "shapefactor" modifier without "constraints", and must survive the +// JSON -> workspace -> JSON round trip with free gammas. The deprecated +// "shapesys" alias is read as well. TEST(HistFactory, HS3ImportShapeFactorModifier) { using namespace RooStats::HistFactory; @@ -745,12 +845,9 @@ TEST(HistFactory, HS3ImportShapeFactorModifier) chan.AddSample(sig); // ShapeFactor on the background: an unconstrained, bin-by-bin scaling. - // Make the gammas constant so that the workspace is well-defined for - // re-export (free shapefactor gammas have no constraints attached). Sample bkg("background", "background", inputFile); ShapeFactor sf; sf.SetName("bkgShape"); - sf.SetConstant(true); bkg.AddShapeFactor(sf); chan.AddSample(bkg); @@ -762,26 +859,30 @@ TEST(HistFactory, HS3ImportShapeFactorModifier) const std::string js = RooJSONFactoryWSTool{*ws}.exportJSONtoString(); - // Rewrite the modifier type for "bkgShape" from "shapesys" to "shapefactor". - // The HistFactory exporter always writes "shapesys", but the importer should - // accept the more accurate "shapefactor" type as well. - const std::regex pattern{"\"name\":\"bkgShape\",\"parameters\":\\[([^\\]]*)\\],\"type\":\"shapesys\""}; - const std::string jsShapeFactor = - std::regex_replace(js, pattern, "\"name\":\"bkgShape\",\"parameters\":[$1],\"type\":\"shapefactor\""); - ASSERT_NE(js, jsShapeFactor) << "Failed to substitute shapesys -> shapefactor in JSON"; + // The free ShapeFactor is written as "shapefactor", without constraints. + const std::regex pattern{"\"name\":\"bkgShape\",\"parameters\":\\[[^\\]]*\\],\"type\":\"shapefactor\""}; + EXPECT_TRUE(std::regex_search(js, pattern)) << js; + EXPECT_EQ(js.find("shapesys"), std::string::npos) << js; RooWorkspace wsFromJson{"new"}; - ASSERT_NO_THROW(RooJSONFactoryWSTool{wsFromJson}.importJSONfromString(jsShapeFactor)) - << "Importer rejected the 'shapefactor' modifier type"; + ASSERT_NO_THROW(RooJSONFactoryWSTool{wsFromJson}.importJSONfromString(js)); - // The imported workspace should expose the same ParamHistFunc gammas. - EXPECT_NE(wsFromJson.var("gamma_bkgShape_bin_0"), nullptr); - EXPECT_NE(wsFromJson.var("gamma_bkgShape_bin_1"), nullptr); + // The imported workspace should expose the same free ParamHistFunc gammas. + ASSERT_NE(wsFromJson.var("gamma_bkgShape_bin_0"), nullptr); + ASSERT_NE(wsFromJson.var("gamma_bkgShape_bin_1"), nullptr); + EXPECT_FALSE(wsFromJson.var("gamma_bkgShape_bin_0")->isConstant()); - // Re-exporting should give back the original JSON, since the writer emits - // type "shapesys" in both cases. + // Re-exporting should give back the original JSON. const std::string js2 = RooJSONFactoryWSTool{wsFromJson}.exportJSONtoString(); EXPECT_EQ(js, js2) << "JSON -> WS -> JSON roundtrip changed the JSON"; + + // The deprecated "shapesys" alias gives the same model. + const std::string jsShapeSys = + std::regex_replace(js, std::regex{"\"type\":\"shapefactor\""}, "\"type\":\"shapesys\""); + ASSERT_NE(js, jsShapeSys); + RooWorkspace wsFromAlias{"alias"}; + ASSERT_NO_THROW(RooJSONFactoryWSTool{wsFromAlias}.importJSONfromString(jsShapeSys)); + EXPECT_EQ(RooJSONFactoryWSTool{wsFromAlias}.exportJSONtoString(), js); } // Issue #20697: Sample::AddShapeFactor() now allows to set the initial value diff --git a/roofit/hs3/inc/RooFitHS3/RooJSONFactoryWSTool.h b/roofit/hs3/inc/RooFitHS3/RooJSONFactoryWSTool.h index c711a16b355cf..ccea81be5c22e 100644 --- a/roofit/hs3/inc/RooFitHS3/RooJSONFactoryWSTool.h +++ b/roofit/hs3/inc/RooFitHS3/RooJSONFactoryWSTool.h @@ -24,6 +24,7 @@ #include #include #include +#include namespace RooFit { namespace JSONIO { @@ -45,6 +46,13 @@ class RooJSONFactoryWSTool { bool allowExportInvalidNames = true; bool allowSanitizeNames = true; bool importNoDomainParametersAsRooConstVars = true; + // Write array distributions (HS3 `array_dist`). If false, no `array_dist` is written: HistFactory constraints + // are named individually, and a RooProdPdf tagged as array distribution is written as `product_dist`. + bool exportArrayDist = true; + // How the constraint terms of HistFactory models enter the likelihood: listed in "aux_distributions", with + // the observed values as constant parameters (false), or listed in "distributions", each paired with a data set + // holding the observed values of its observables (true). + bool exportConstraintData = false; }; static Config &config(); @@ -200,6 +208,14 @@ class RooJSONFactoryWSTool { } std::string exportTransformed(const RooAbsReal *original, const std::string &suffix, const std::string &formula); + std::string exportAsArrayComponent(RooAbsPdf const &pdf, std::string const &arrayName); + /// The name of the top-level object whose export is in progress (the root of the current export recursion). + std::string const &exportRootName() const { return _exportRootName; } + void recordConstraintTerms(std::string const &pdfName, std::vector const &terms); + /// Record that `observable` is the observable of the constraint term `pdf` (the side that is not the constrained + /// parameter). A Gaussian whose mean it is gets written with x and mean swapped, so that the observable is x. + void setConstraintObservable(RooAbsPdf const &pdf, std::string const &observable); + std::string const *constraintObservable(std::string const &pdfName) const; void setAttribute(const std::string &obj, const std::string &attrib); bool hasAttribute(const std::string &obj, const std::string &attrib); @@ -212,11 +228,11 @@ class RooJSONFactoryWSTool { private: template T *requestImpl(const std::string &objname); - void exportObject(RooAbsArg const &func, std::set &exportedObjectNames); + void exportObject(RooAbsArg const &func); // To export multiple objects sorted alphabetically template - void exportObjects(T const &args, std::set &exportedObjectNames) + void exportObjects(T const &args) { RooArgSet argSet; for (RooAbsArg const *arg : args) { @@ -224,9 +240,12 @@ class RooJSONFactoryWSTool { } argSet.sort(); for (RooAbsArg *arg : argSet) { - exportObject(*arg, exportedObjectNames); + exportObject(*arg); } } + void exportArrayDists(); + bool exportConstraintData(std::string const &term, RooArgSet const *globalObservables, + RooFit::Detail::JSONNode &rootnode); void exportData(RooAbsData const &data); RooJSONFactoryWSTool::CombinedData exportCombinedData(RooAbsData const &data); @@ -261,6 +280,17 @@ class RooJSONFactoryWSTool { std::unique_ptr _domains; std::vector _serversToExport; std::vector _serversToDelete; + // Names of the objects written so far, to export each object once. + std::set _exportedObjectNames; + // Array distributions to write at the end of the export: their components, and the array of each component. + std::map> _arrayComponents; + std::map _arrayOfComponent; + std::string _exportRootName; + int _exportDepth = 0; + // The constraint terms of the exported HistFactory channels, as written, and the data sets written for them. + std::map> _constraintTermsOfPdf; + std::set _writtenConstraintData; + std::map _constraintObservables; // Name-keyed indices over the top-level "functions" and "distributions" // sequences of the input JSON. Built once at the start of importAllNodes() diff --git a/roofit/hs3/src/Domains.cxx b/roofit/hs3/src/Domains.cxx index a68513335d094..3894ad52ae32a 100644 --- a/roofit/hs3/src/Domains.cxx +++ b/roofit/hs3/src/Domains.cxx @@ -15,6 +15,7 @@ #include #include +#include #include #include #include @@ -22,6 +23,8 @@ #include +#include + namespace RooFit { namespace JSONIO { namespace Detail { @@ -30,17 +33,63 @@ constexpr static auto defaultDomainName = "default_domain"; namespace { -void writeBound(RooFit::Detail::JSONNode &node, double value) +// JSON can't represent infinity, so an infinite bound is expressed by omitting it. An absent bound is read as +// infinite by readBound() in JSONIOUtils. +void writeBound(RooFit::Detail::JSONNode &node, const char *key, double value) { - if (RooNumber::isInfinite(value)) { - node.set_null(); - } else { - node << value; + if (!RooNumber::isInfinite(value)) { + node[key] << value; + } +} + +// Read the "categories" of a categorical axis. The labels must be unique non-empty strings, and the axis must not have +// the fields of a continuous axis. +std::vector readCategories(RooFit::Detail::JSONNode const &axis, std::string const &context) +{ + for (const char *key : {"min", "max", "nbins", "edges"}) { + if (axis.has_child(key)) { + RooJSONFactoryWSTool::error(context + " is categorical and must not have \"" + key + "\""); + } + } + auto const &node = axis["categories"]; + if (!node.is_seq() || node.num_children() == 0) { + RooJSONFactoryWSTool::error(context + ": \"categories\" must be a non-empty list of strings"); + } + std::vector labels; + std::set seen; + for (auto const &label : node.children()) { + if (label.is_seq() || label.is_map() || label.is_number() || label.val().empty()) { + RooJSONFactoryWSTool::error(context + ": the entries of \"categories\" must be non-empty strings"); + } + if (!seen.insert(label.val()).second) { + RooJSONFactoryWSTool::error(context + ": the category \"" + label.val() + "\" is listed twice"); + } + labels.push_back(label.val()); } + return labels; } } // namespace +std::vector const *Domains::categories(const char *name) const +{ + for (auto const &domain : _map) { + if (auto const *labels = domain.second.categories(name)) { + return labels; + } + } + return nullptr; +} + +std::vector const *Domains::ProductDomain::categories(const char *name) const +{ + auto found = _map.find(name); + if (found == _map.end() || found->second.categories.empty()) { + return nullptr; + } + return &found->second.categories; +} + void Domains::populate(RooWorkspace &ws) const { auto default_domain = _map.find(defaultDomainName); @@ -186,6 +235,10 @@ void Domains::ProductDomain::writeVariable(RooRealVar &var) const auto found = _map.find(var.GetName()); if (found != _map.end()) { auto const &elem = found->second; + if (!elem.categories.empty()) { + RooJSONFactoryWSTool::error(std::string("the variable '") + var.GetName() + + "' is continuous, but a domain lists it as categorical"); + } if (elem.hasMin) { if (RooNumber::isInfinite(elem.min)) { var.removeMin(); @@ -215,7 +268,22 @@ void Domains::ProductDomain::readJSON(RooFit::Detail::JSONNode const &node) RooJSONFactoryWSTool::error("only domains of type \"product_domain\" are currently supported!"); } for (auto const &varNode : node["axes"].children()) { - auto &elem = _map[RooJSONFactoryWSTool::name(varNode)]; + const std::string varName = RooJSONFactoryWSTool::name(varNode); + auto &elem = _map[varName]; + + if (varNode.has_child("categories")) { + std::vector labels = readCategories(varNode, "the axis '" + varName + "' of a domain"); + if (!elem.categories.empty() && elem.categories != labels) { + RooJSONFactoryWSTool::error("the categorical axis '" + varName + + "' is listed twice in a domain, with different categories"); + } + elem.categories = std::move(labels); + continue; + } + if (!elem.categories.empty()) { + RooJSONFactoryWSTool::error("the axis '" + varName + + "' is listed in a domain as categorical and as continuous"); + } if (varNode.has_child("min")) { elem.min = readBound(varNode, "min", -RooNumber::infinity()); @@ -259,17 +327,30 @@ void Domains::ProductDomain::writeJSON(RooFit::Detail::JSONNode &node) const for (auto const &item : _map) { auto const &elem = item.second; RooFit::Detail::JSONNode &varnode = RooJSONFactoryWSTool::appendNamedChild(variablesNode, item.first); - writeBound(varnode["min"], elem.hasMin ? elem.min : -RooNumber::infinity()); - writeBound(varnode["max"], elem.hasMax ? elem.max : RooNumber::infinity()); - writeBinning(varnode, elem); + if (!elem.categories.empty()) { + varnode["categories"].fill_seq(elem.categories); + continue; + } + const double min = elem.hasMin ? elem.min : -RooNumber::infinity(); + const double max = elem.hasMax ? elem.max : RooNumber::infinity(); + writeBound(varnode, "min", min); + writeBound(varnode, "max", max); + // A binning requires a bounded range. + if (!RooNumber::isInfinite(min) && !RooNumber::isInfinite(max)) { + writeBinning(varnode, elem); + } } } void Domains::ProductDomain::populate(RooWorkspace &ws) const { for (auto const &item : _map) { const auto &name = item.first; + const auto &elem = item.second; + if (!elem.categories.empty()) { + populateCategory(name, elem.categories, ws); + continue; + } if (!ws.arg(name)) { - const auto &elem = item.second; const double vMin = elem.hasMin ? elem.min : -RooNumber::infinity(); const double vMax = elem.hasMax ? elem.max : RooNumber::infinity(); RooRealVar var{name.c_str(), name.c_str(), vMin, vMax}; @@ -278,9 +359,40 @@ void Domains::ProductDomain::populate(RooWorkspace &ws) const } } } +// Create a RooCategory for a categorical axis, with one state per category. The states are defined in the order of the +// document, with the position as index, so that the bins of a RooDataHist follow the order of "categories". An +// existing RooCategory of that name must have the same labels. +void Domains::ProductDomain::populateCategory(std::string const &name, std::vector const &labels, + RooWorkspace &ws) +{ + if (RooAbsArg *existing = ws.arg(name)) { + auto *cat = dynamic_cast(existing); + if (!cat) { + RooJSONFactoryWSTool::error("the axis '" + name + "' is categorical, but '" + name + + "' already exists and is not a category"); + } + std::set existingLabels; + for (auto const &state : *cat) { + existingLabels.insert(state.first); + } + if (existingLabels != std::set(labels.begin(), labels.end())) { + RooJSONFactoryWSTool::error("the categories of the axis '" + name + + "' differ from the states of the existing category '" + name + "'"); + } + return; + } + RooCategory cat{name.c_str(), name.c_str()}; + for (std::size_t i = 0; i < labels.size(); ++i) { + cat.defineType(labels[i], static_cast(i)); + } + ws.import(cat, RooFit::Silence()); +} + void Domains::ProductDomain::registerBinnings(const char *name, RooWorkspace &ws) const { for (auto const &item : _map) { + if (!item.second.categories.empty()) + continue; auto *var = ws.var(item.first); if (!var) continue; diff --git a/roofit/hs3/src/Domains.h b/roofit/hs3/src/Domains.h index 1c39a52ffa206..a307ae5ac8009 100644 --- a/roofit/hs3/src/Domains.h +++ b/roofit/hs3/src/Domains.h @@ -41,6 +41,7 @@ class Domains { void writeJSON(RooFit::Detail::JSONNode &) const; bool hasVariable(const char *name) const; + std::vector const *categories(const char *name) const; void populate(RooWorkspace &ws) const; @@ -54,6 +55,7 @@ class Domains { void writeJSON(RooFit::Detail::JSONNode &) const; bool hasVariable(const char *name) const; + std::vector const *categories(const char *name) const; void populate(RooWorkspace &ws) const; void registerBinnings(const char *name, RooWorkspace &ws) const; @@ -67,8 +69,11 @@ class Domains { bool hasNBins = false; int nBins = 0; std::vector edges; + // The labels of a categorical axis, in the order given in the document. Empty for continuous axes. + std::vector categories; }; + static void populateCategory(std::string const &name, std::vector const &labels, RooWorkspace &ws); static void applyBinning(RooRealVar &var, ProductDomainElement const &elem, const char *name = nullptr); static void readBinning(ProductDomainElement &elem, RooAbsBinning const &binning); static void writeBinning(RooFit::Detail::JSONNode &node, ProductDomainElement const &elem); diff --git a/roofit/hs3/src/JSONFactories_HistFactory.cxx b/roofit/hs3/src/JSONFactories_HistFactory.cxx index c5f34761b65f7..c3f6bb4d72337 100644 --- a/roofit/hs3/src/JSONFactories_HistFactory.cxx +++ b/roofit/hs3/src/JSONFactories_HistFactory.cxx @@ -444,6 +444,55 @@ void exportGammaParameters(JSONNode &mod, std::vector const ¶m mod["parameters"].fill_seq(paramNames); } +// The names under which the constraint distributions of a model are written: the array distribution that holds +// the constraint, or its own name. +using ConstraintNames = std::map; + +/// Export the names of the constraint distributions of a modifier, with null for unconstrained parameters. Nothing is +/// written if all parameters are unconstrained. +void exportConstraintNames(JSONNode &mod, std::vector const &constraints, + ConstraintNames const &names) +{ + if (std::none_of(constraints.begin(), constraints.end(), [](auto *pdf) { return pdf != nullptr; })) { + return; + } + auto &out = mod["constraints"].set_seq(); + for (auto *constraint : constraints) { + if (!constraint) { + out.append_child().set_null(); + continue; + } + auto found = names.find(constraint); + out.append_child() << (found != names.end() ? found->second : std::string(constraint->GetName())); + } +} + +// Register the constraint `constraint` of a modifier parameter for the export and return the name under which it is +// written. Constraints are collected into one array distribution per type and top-level pdf (or, if the constraint +// already belongs to an array distribution, that one), unless the arrays are switched off or the constraint cannot be +// a component. In those cases it is written on its own. +std::string exportConstraint(RooJSONFactoryWSTool &tool, RooAbsPdf const &constraint) +{ + if (RooJSONFactoryWSTool::config().exportArrayDist) { + if (auto info = arrayComponentInfo(constraint)) { + std::string arrayName; + for (RooAbsArg const *client : constraint.clients()) { + if (client->getAttribute(arrayDistTag)) { + arrayName = client->GetName(); + break; + } + } + if (arrayName.empty()) { + // "poisson_dist" -> "_poisson_constraints" + arrayName = tool.exportRootName() + "_" + removeSuffix(info->type, "_dist") + "_constraints"; + } + return tool.exportAsArrayComponent(constraint, arrayName); + } + } + tool.queueExport(constraint); + return constraint.GetName(); +} + RooRealVar &createNominal(RooWorkspace &ws, std::string const &parname, double val, double min, double max) { RooRealVar &nom = getOrCreate(ws, "nom_" + parname, val, min, max); @@ -481,13 +530,63 @@ RooAbsPdf *findNamedConstraint(RooJSONFactoryWSTool &tool, std::string const &co return nullptr; } +// A description of a modifier for error messages. +std::string modifierContext(const JSONNode &mod, std::string const &sample) +{ + std::string out = "the " + mod["type"].val() + " modifier"; + if (mod.has_child("name")) { + out += " '" + mod["name"].val() + "'"; + } + return out + " of sample '" + sample + "'"; +} + +// The component of the array distribution `array` that constrains `param`: the one that depends on it. +RooAbsPdf &arrayComponentFor(RooProdPdf const &array, RooAbsArg const ¶m, std::string const &context) +{ + std::vector matches; + for (RooAbsPdf *pdf : static_range_cast(array.pdfList())) { + if (pdf->dependsOn(param)) { + matches.push_back(pdf); + } + } + if (matches.size() == 1) { + return *matches.front(); + } + const std::string what = "the array distribution '" + std::string(array.GetName()) + + "' named as the constraint of '" + param.GetName() + "' (" + context + ") has "; + if (matches.empty()) { + RooJSONFactoryWSTool::error(what + "no component that depends on this parameter"); + } + RooJSONFactoryWSTool::error(what + std::to_string(matches.size()) + + " components that depend on this parameter, but a parameter can only be constrained " + "by one component"); +} + +// Find the constraint distribution `constraintName` of the parameter `param`. If the name resolves to an array +// distribution, the component that depends on the parameter is returned: that is the term that enters the likelihood +// for this parameter, and the whole array enters once through the components of all of its parameters. Returns +// nullptr if the name cannot be resolved. +RooAbsPdf *findNamedConstraint(RooJSONFactoryWSTool &tool, std::string const &constraintName, std::string const &sample, + RooAbsArg const ¶m, std::string const &context) +{ + RooAbsPdf *constraint = findNamedConstraint(tool, constraintName, sample); + if (!constraint || !constraint->getAttribute(arrayDistTag)) { + return constraint; + } + auto const *array = dynamic_cast(constraint); + if (!array) { + return constraint; + } + return &arrayComponentFor(*array, param, context); +} + RooAbsPdf &createLegacyConstraint(RooJSONFactoryWSTool &tool, const JSONNode &mod, RooRealVar ¶m, std::string const &constraintType) { if (constraintType == "Gauss") { param.setError(1.0); - return getOrCreate(*tool.workspace(), constraintName(param.GetName()), param, - *tool.workspace()->var(std::string("nom_") + param.GetName()), 1.); + RooRealVar &nominal = createNominal(*tool.workspace(), param.GetName(), 0.0, -10, 10); + return getOrCreate(*tool.workspace(), constraintName(param.GetName()), param, nominal, 1.); } RooJSONFactoryWSTool::error("legacy constraint value '" + constraintType + "' for modifier '" + @@ -514,7 +613,10 @@ ParamHistFunc &createPHF(const std::string &phfname, std::string const &sysname, gammas.add(getOrCreate(ws, name, 1., gammaMin, gammaMax)); } - auto &phf = tool.wsEmplace(phfname, observables, gammas); + // The gammas and `vals` are given in the order of HS3 (last axis fastest), which is also the order in which they + // are paired below. ParamHistFunc expects its parameters in TH1 order. + checkParamHistFuncOrderIsDefined(tool.hs3VersionInput(), observables, "the modifier '" + sysname + "'"); + auto &phf = tool.wsEmplace(phfname, observables, toParamHistFuncOrder(gammas, observables)); if (vals.size() > 0) { if (!createConstraints) { @@ -553,13 +655,77 @@ const JSONNode *findStaterror(const JSONNode &comp) return nullptr; } +// The entries of the HS3 "constraints" of a modifier: the names of the constraint distributions, or std::nullopt for +// unconstrained parameters (null entries). +std::vector> constraintNames(const JSONNode &mod) +{ + std::vector> out; + auto *constraints = mod.find("constraints"); + if (!constraints) { + return out; + } + if (!constraints->is_seq()) { + RooJSONFactoryWSTool::error("\"constraints\" of modifier '" + RooJSONFactoryWSTool::name(mod) + + "' must be a list"); + } + for (const auto &entry : constraints->children()) { + out.push_back(entry.is_null() ? std::nullopt : std::make_optional(entry.val())); + } + return out; +} + +// The parameter of a modifier with a single parameter, given in the HS3 "parameters" or in the legacy "parameter". +std::optional modifierParameter(const JSONNode &mod) +{ + if (auto *parameters = mod.find("parameters")) { + if (!parameters->is_seq() || parameters->num_children() != 1) { + RooJSONFactoryWSTool::error("\"parameters\" of the " + mod["type"].val() + + " modifier must be a list with one entry"); + } + return parameters->child(0).val(); + } + if (auto *parameter = mod.find("parameter")) { + return parameter->val(); + } + return std::nullopt; +} + +// The constraint distribution named in the HS3 "constraints" of a modifier with a single parameter. Returns an empty +// string if the parameter is explicitly unconstrained, and std::nullopt if the modifier has no "constraints". +std::optional singleConstraintName(const JSONNode &mod) +{ + if (!mod.has_child("constraints")) { + return std::nullopt; + } + const auto names = constraintNames(mod); + if (names.size() != 1) { + RooJSONFactoryWSTool::error("\"constraints\" of the " + mod["type"].val() + + " modifier must have one entry per parameter"); + } + return names.front().value_or(""); +} + RooAbsPdf & getOrCreateConstraint(RooJSONFactoryWSTool &tool, const JSONNode &mod, RooRealVar ¶m, const std::string &sample) { + const std::string context = modifierContext(mod, sample); + // HS3 references the constraint distributions by name. + if (auto name = singleConstraintName(mod); name && !name->empty()) { + auto *constraint = findNamedConstraint(tool, *name, sample, param, context); + if (!constraint) { + RooJSONFactoryWSTool::error("unable to find the constraint distribution '" + *name + "' of the modifier of '" + + param.GetName() + "'"); + } + if (auto gauss = dynamic_cast(constraint)) { + param.setError(gauss->getSigma().getVal()); + } + return *constraint; + } + JSONNode const *constrName = mod.find("constraint_name"); if (constrName) { auto constraint_name = constrName->val(); - auto constraint = findNamedConstraint(tool, constraint_name, sample); + auto constraint = findNamedConstraint(tool, constraint_name, sample, param, context); if (!constraint) { RooJSONFactoryWSTool::error("unable to find definition of of constraint '" + constraint_name + "' for modifier '" + RooJSONFactoryWSTool::name(mod) + "'"); @@ -572,7 +738,7 @@ getOrCreateConstraint(RooJSONFactoryWSTool &tool, const JSONNode &mod, RooRealVa if (auto constr = mod.find("constraint")) { std::string constraintValue = constr->val(); - if (auto *constraint = findNamedConstraint(tool, constraintValue, sample)) { + if (auto *constraint = findNamedConstraint(tool, constraintValue, sample, param, context)) { if (auto gauss = dynamic_cast(constraint)) { param.setError(gauss->getSigma().getVal()); } @@ -600,6 +766,9 @@ getOrCreateConstraint(RooJSONFactoryWSTool &tool, const JSONNode &mod, RooRealVa } bool hasConstraint(const JSONNode &mod) { + if (auto name = singleConstraintName(mod)) { + return !name->empty(); + } return mod.has_child("constraint") || mod.has_child("constraint_name") || mod.has_child("constraint_type"); } @@ -618,7 +787,6 @@ void addSysConstraint(RooJSONFactoryWSTool &tool, const JSONNode &mod, RooRealVa } return; } - createNominal(*tool.workspace(), param.GetName(), 0.0, -10, 10); constraints.add(getOrCreateConstraint(tool, mod, param, sample)); } @@ -671,7 +839,7 @@ bool importHistSample(RooJSONFactoryWSTool &tool, RooDataHist &dh, RooArgSet con RooAbsArg const *mcStatObject, const std::string &fprefix, const JSONNode &p, const std::optional &defaultInterpolation, LegacyInterpolationDefault &legacyInterpolationDefault, bool &warnedUnconstrained, - RooArgSet &constraints) + bool &warnedShapeSys, RooArgSet &constraints) { RooWorkspace &ws = *tool.workspace(); @@ -711,17 +879,17 @@ bool importHistSample(RooJSONFactoryWSTool &tool, RooDataHist &dh, RooArgSet con int idx = 0; for (const auto &mod : p["modifiers"].children()) { std::string const &modtype = mod["type"].val(); + const std::optional parameter = + modtype == "normfactor" || modtype == "normsys" || modtype == "histosys" ? modifierParameter(mod) + : std::nullopt; std::string const &sysname = - mod.has_child("name") - ? mod["name"].val() - : (mod.has_child("parameter") ? mod["parameter"].val() : "syst_" + std::to_string(idx)); + mod.has_child("name") ? mod["name"].val() : parameter.value_or("syst_" + std::to_string(idx)); ++idx; if (modtype == "staterror") { // this is dealt with at a different place, ignore it for now } else if (modtype == "normfactor") { - // The "name" is only a label of the modifier, the parameter is given by "parameter" if present. - auto *parameter = mod.find("parameter"); - RooRealVar &constrParam = getOrCreate(ws, parameter ? parameter->val() : sysname, 1., -3, 5); + // The "name" is only a label of the modifier, the parameter is given by "parameters" if present. + RooRealVar &constrParam = getOrCreate(ws, parameter.value_or(sysname), 1., -3, 5); constrParam.setError(0.0); normElems.add(constrParam); if (hasConstraint(mod)) { @@ -729,8 +897,7 @@ bool importHistSample(RooJSONFactoryWSTool &tool, RooDataHist &dh, RooArgSet con constraints.add(getOrCreateConstraint(tool, mod, constrParam, sampleName)); } } else if (modtype == "normsys") { - auto *parameter = mod.find("parameter"); - std::string parname(parameter ? parameter->val() : "alpha_" + sysname); + std::string parname(parameter.value_or("alpha_" + sysname)); auto &par = getOrCreate(ws, parname, 0., -5, 5); overall_nps.add(par); auto &data = mod["data"]; @@ -756,8 +923,7 @@ bool importHistSample(RooJSONFactoryWSTool &tool, RooDataHist &dh, RooArgSet con addSysConstraint(tool, mod, par, sysname, sampleName, warnedUnconstrained, constraints); } else if (modtype == "histosys") { - auto *parameter = mod.find("parameter"); - std::string parname(parameter ? parameter->val() : "alpha_" + sysname); + std::string parname(parameter.value_or("alpha_" + sysname)); auto &par = getOrCreate(ws, parname, 0., -5, 5); histNps.add(par); auto &data = mod["data"]; @@ -773,8 +939,13 @@ bool importHistSample(RooJSONFactoryWSTool &tool, RooDataHist &dh, RooArgSet con legacyInterpolationDefault)); addSysConstraint(tool, mod, par, sysname, sampleName, warnedUnconstrained, constraints); } else if (modtype == "shapesys" || modtype == "shapefactor") { - std::string funcName = channelName + "_" + sysname + "_ShapeSys"; - // funcName should be "__ShapeSys" + // HS3 has a single per-bin modifier type, "shapefactor", whose parameters are constrained by the + // distributions in "constraints", if given. "shapesys" is a deprecated alias. + if (modtype == "shapesys" && !isHS3VersionBefore(tool.hs3VersionInput(), 0, 3) && !warnedShapeSys) { + RooJSONFactoryWSTool::warning("RooFitHS3: the modifier type \"shapesys\" is deprecated, please use " + "\"shapefactor\" instead."); + warnedShapeSys = true; + } std::vector vals; if (mod["data"].has_child("vals")) { for (const auto &v : mod["data"]["vals"].children()) { @@ -790,26 +961,12 @@ bool importHistSample(RooJSONFactoryWSTool &tool, RooDataHist &dh, RooArgSet con "' with neither values nor parameters!"); } std::string constraint = "unknown"; - std::vector constraintPdfs; bool const hasConstraintList = mod.has_child("constraints"); if (hasConstraintList) { - for (const auto &v : mod["constraints"].children()) { - if (v.is_null()) { - constraintPdfs.push_back(nullptr); - } else { - std::string constraintName = v.val(); - auto *constraintPdf = findNamedConstraint(tool, constraintName, sampleName); - if (!constraintPdf) { - RooJSONFactoryWSTool::error("unable to find definition of constraint '" + constraintName + - "' for modifier '" + RooJSONFactoryWSTool::name(mod) + "'"); - } - constraintPdfs.push_back(constraintPdf); - } - } std::size_t const nGammas = std::max(vals.size(), parnames.size()); - if (constraintPdfs.size() != nGammas) { + if (mod["constraints"].num_children() != nGammas) { std::stringstream ss; - ss << "modifier '" << RooJSONFactoryWSTool::name(mod) << "' has " << constraintPdfs.size() + ss << "modifier '" << RooJSONFactoryWSTool::name(mod) << "' has " << mod["constraints"].num_children() << " constraints, but " << nGammas << " parameters"; RooJSONFactoryWSTool::error(ss.str()); } @@ -826,11 +983,48 @@ bool importHistSample(RooJSONFactoryWSTool &tool, RooDataHist &dh, RooArgSet con "not a supported legacy constraint type"); } } - shapeElems.add(createPHF(funcName, sysname, parnames, vals, tool, constraints, varlist, constraint, - defaultGammaMin, defaultShapeSysGammaMax, minShapeUncertainty, - /*createConstraints=*/!hasConstraintList)); - for (auto *constraintPdf : constraintPdfs) { - if (constraintPdf) { + // Without "constraints", only the legacy constraint keys, or documents before HS3 0.3, make RooFit + // create constraints from "data.vals". Otherwise, the parameters are free. + const bool legacyConstraint = mod.has_child("constraint_type") || mod.has_child("constraint"); + const bool createConstraints = + !hasConstraintList && (legacyConstraint || isHS3VersionBefore(tool.hs3VersionInput(), 0, 3)); + if (!hasConstraintList && !createConstraints && !vals.empty()) { + oocoutW(nullptr, ObjectHandling) + << "RooFitHS3: the \"data.vals\" of modifier '" << sysname << "' of sample '" << sampleName + << "' are ignored, because it names no constraints, so its " + << "parameters are free" << std::endl; + vals.clear(); + } + bool constrained = createConstraints && !vals.empty() && constraint != "Const"; + if (hasConstraintList) { + for (const auto &v : mod["constraints"].children()) { + constrained |= !v.is_null(); + } + } + // Follow the naming of HistFactory: "__ShapeSys" for constrained and + // "__shapeFactor" for free parameters. + const std::string funcName = channelName + "_" + sysname + (constrained ? "_ShapeSys" : "_shapeFactor"); + ParamHistFunc &phf = + createPHF(funcName, sysname, parnames, vals, tool, constraints, varlist, constraint, defaultGammaMin, + defaultShapeSysGammaMax, minShapeUncertainty, createConstraints); + shapeElems.add(phf); + if (hasConstraintList) { + // The constraints are resolved per gamma (in bin order, like "constraints"), which also picks the + // right component out of an array distribution. + const std::vector gammas = paramHistFuncParamsInBinOrder(phf); + std::size_t i = 0; + for (const auto &v : mod["constraints"].children()) { + RooAbsReal const &gamma = *gammas[i++]; + if (v.is_null()) { + continue; + } + std::string constraintName = v.val(); + auto *constraintPdf = + findNamedConstraint(tool, constraintName, sampleName, gamma, modifierContext(mod, sampleName)); + if (!constraintPdf) { + RooJSONFactoryWSTool::error("unable to find definition of constraint '" + constraintName + + "' for modifier '" + RooJSONFactoryWSTool::name(mod) + "'"); + } constraints.add(*constraintPdf); } } @@ -901,6 +1095,13 @@ class HistFactoryImporter : public RooFit::JSONIO::Importer { std::vector sumW2; std::vector gammaParnames; std::string modifierStatErrType; + std::optional>> statErrConstraintNames; + for (auto const &axis : p["axes"].children()) { + if (axis.has_child("categories")) { + RooJSONFactoryWSTool::error("the axis '" + axis["name"].val() + "' of '" + name + + "' is categorical, but the axes of a HistFactory prediction are binned"); + } + } RooArgSet observables = RooJSONFactoryWSTool::readAxes(p); std::string fprefix = name; @@ -916,9 +1117,13 @@ class HistFactoryImporter : public RooFit::JSONIO::Importer { sumW.resize(nbins); sumW2.resize(nbins); } + // A sample without "errors" contributes no statistical uncertainty in HS3. Before HS3 0.3, RooFit + // assumed Poisson uncertainties for it. + const bool useWeightSquared = + comp["data"].has_child("errors") || isHS3VersionBefore(tool->hs3VersionInput(), 0, 3); for (size_t i = 0; i < nbins; ++i) { sumW[i] += dh->weight(i); - sumW2[i] += dh->weightSquared(i); + sumW2[i] += useWeightSquared ? dh->weightSquared(i) : 0.; } if (gammaParnames.empty()) { if (auto staterrorParams = staterror->find("parameters")) { @@ -927,7 +1132,12 @@ class HistFactoryImporter : public RooFit::JSONIO::Importer { } } } - // The constraint type given on the staterror modifier takes precedence over the channel-level setting. + // HS3 names the constraint distributions of the gammas in "constraints". + if (!statErrConstraintNames && staterror->has_child("constraints")) { + statErrConstraintNames = constraintNames(*staterror); + } + // The legacy constraint type given on the staterror modifier takes precedence over the channel-level + // setting. if (modifierStatErrType.empty()) { // ROOT 6.38 and 6.40.00 wrote "constraint_type". auto constraint = staterror->find("constraint"); @@ -950,6 +1160,11 @@ class HistFactoryImporter : public RooFit::JSONIO::Importer { if (!modifierStatErrType.empty()) { statErrType = modifierStatErrType; } + // Without any constraint information, the gammas are unconstrained in HS3. Older documents relied on RooFit + // creating Poisson constraints. + const bool legacyStatErrConstraints = statErrConstraintNames || !modifierStatErrType.empty() || + p.has_child(::Literals::staterror) || + isHS3VersionBefore(tool->hs3VersionInput(), 0, 3); RooAbsArg *mcStatObject = nullptr; RooArgSet constraints; @@ -968,21 +1183,45 @@ class HistFactoryImporter : public RooFit::JSONIO::Importer { errs[i] = std::max(errs[i], 0.); } - mcStatObject = - &createPHF("mc_stat_" + channelName, "stat_" + channelName, gammaParnames, errs, *tool, constraints, - observables, statErrType, defaultGammaMin, defaultStatErrorGammaMax, statErrThresh); + const bool createConstraints = !statErrConstraintNames && legacyStatErrConstraints; + mcStatObject = &createPHF("mc_stat_" + channelName, "stat_" + channelName, gammaParnames, errs, *tool, + constraints, observables, statErrType, defaultGammaMin, defaultStatErrorGammaMax, + statErrThresh, createConstraints); + if (statErrConstraintNames) { + if (statErrConstraintNames->size() != errs.size()) { + RooJSONFactoryWSTool::error("the staterror modifier of channel '" + name + "' has " + + std::to_string(statErrConstraintNames->size()) + " constraints, but " + + std::to_string(errs.size()) + " parameters"); + } + const std::vector gammas = + paramHistFuncParamsInBinOrder(static_cast(*mcStatObject)); + for (std::size_t i = 0; i < statErrConstraintNames->size(); ++i) { + auto const &constraintName = (*statErrConstraintNames)[i]; + if (!constraintName) { + continue; + } + auto *constraint = findNamedConstraint(*tool, *constraintName, name, *gammas[i], + "the staterror modifier of channel '" + name + "'"); + if (!constraint) { + RooJSONFactoryWSTool::error("unable to find the constraint distribution '" + *constraintName + + "' of the staterror modifier of channel '" + name + "'"); + } + constraints.add(*constraint); + } + } } LegacyInterpolationDefault legacyInterpolationDefault; legacyInterpolationDefault.allowed = isHS3VersionBefore(tool->hs3VersionInput(), 0, 3); bool warnedUnconstrained = false; + bool warnedShapeSys = false; int idx = 0; RooArgList funcs; RooArgList coefs; for (const auto &comp : p["samples"].children()) { importHistSample(*tool, *data[idx], observables, mcStatObject, fprefix, comp, defaultInterpolation, - legacyInterpolationDefault, warnedUnconstrained, constraints); + legacyInterpolationDefault, warnedUnconstrained, warnedShapeSys, constraints); ++idx; std::string const &compName = RooJSONFactoryWSTool::name(comp); @@ -1160,9 +1399,9 @@ struct HistoSys { high.assign(h->dataHist().weightArray(), h->dataHist().weightArray() + h->dataHist().numEntries()); } }; +// A per-bin modifier ("shapefactor"): HistFactory's ShapeSys (constrained) or ShapeFactor (free). struct ShapeSys { std::string name; - std::vector constraints; std::vector constraintPdfs; std::vector parameters; ShapeSys(const std::string &n) : name{n} {} @@ -1323,23 +1562,18 @@ struct Sample { std::vector otherElements; bool useBarlowBeestonLight = false; std::vector staterrorParameters; + std::vector staterrorConstraints; Sample(const std::string &n) : name{n} {} }; -void addNormFactor(RooRealVar const *par, Sample &sample, RooWorkspace *ws) +void addNormFactor(RooRealVar const *par, Sample &sample) { - std::string parname = par->GetName(); - bool isConstrained = false; - for (RooAbsArg const *pdf : ws->allPdfs()) { - if (auto gauss = dynamic_cast(pdf)) { - if (parname == gauss->getX().GetName()) { - sample.normfactors.emplace_back(*par, gauss); - isConstrained = true; - } - } - } - if (!isConstrained) + // The constraint is found among the clients of the parameter, whichever of x and mean of a Gaussian it is. + if (RooAbsPdf const *constraint = findConstraint(const_cast(par))) { + sample.normfactors.emplace_back(*par, constraint); + } else { sample.normfactors.emplace_back(*par); + } } struct Channel { @@ -1348,7 +1582,6 @@ struct Channel { std::map tot_yield; std::map tot_yield2; std::map rel_errors; - std::string statErrorConstraintType; RooArgSet const *varSet = nullptr; long unsigned int nBins = 0; }; @@ -1378,6 +1611,19 @@ Channel readChannel(RooJSONFactoryWSTool *tool, const std::string &pdfname, cons if (channel.varSet == nullptr) { channel.varSet = dataHist.get(); channel.nBins = dataHist.numEntries(); + } else { + // The bin contents of all samples are written in the order of the axes of the channel (the last axis + // varies fastest), which is the order of the observables of the first histogram. + RooArgSet const &vars = *dataHist.get(); + bool sameOrder = vars.size() == channel.varSet->size(); + for (std::size_t i = 0; sameOrder && i < vars.size(); ++i) { + sameOrder = std::string(vars[i]->GetName()) == (*channel.varSet)[i]->GetName(); + } + if (!sameOrder) { + RooJSONFactoryWSTool::error("the histograms of channel '" + channel.name + + "' don't have their observables in the same order, so their bin " + "contents can't be written along the same axes"); + } } if (sample.hist.empty()) { auto *w = dataHist.weightArray(); @@ -1399,7 +1645,7 @@ Channel readChannel(RooJSONFactoryWSTool *tool, const std::string &pdfname, cons sample.normfactors.emplace_back(*constVar); } } else if (auto par = dynamic_cast(e)) { - addNormFactor(par, sample, ws); + addNormFactor(par, sample); } else if (auto hf = dynamic_cast(e)) { updateObservables(hf->dataHist()); } else if (ParamHistFunc *phf = dynamic_cast(e); phf && allRooRealVar(phf->paramList())) { @@ -1508,44 +1754,35 @@ Channel readChannel(RooJSONFactoryWSTool *tool, const std::string &pdfname, cons for (ParamHistFunc *phf : phfs) { if (startsWith(std::string(phf->GetName()), "mc_stat_")) { // MC stat uncertainty int idx = 0; - for (const auto &g : phf->paramList()) { + // In the order of the bins (last axis fastest), like the contents of the samples. + for (RooAbsReal *g : paramHistFuncParamsInBinOrder(*phf)) { sample.staterrorParameters.push_back(static_cast(g)); ++idx; RooAbsPdf *constraint = findConstraint(g); + sample.staterrorConstraints.push_back(constraint); channel.tot_yield[idx] += sample.hist[idx - 1]; channel.tot_yield2[idx] += (sample.hist[idx - 1] * sample.hist[idx - 1]); if (constraint) { channel.rel_errors[idx] = constraintRelError(*constraint, *g); - channel.statErrorConstraintType = dynamic_cast(constraint) ? "Poisson" : "Gauss"; } } sample.useBarlowBeestonLight = true; - } else { // other ShapeSys + } else { // ShapeSys or ShapeFactor ShapeSys sys(phf->GetName()); erasePrefix(sys.name, channel.name + "_"); - bool isshapesys = eraseSuffix(sys.name, "_ShapeSys") || eraseSuffix(sys.name, "_shapeSys"); - bool isshapefactor = eraseSuffix(sys.name, "_ShapeFactor") || eraseSuffix(sys.name, "_shapeFactor"); + for (std::string_view suffix : {"_ShapeSys", "_shapeSys", "_ShapeFactor", "_shapeFactor"}) { + if (eraseSuffix(sys.name, suffix)) + break; + } - for (const auto &g : phf->paramList()) { + // Both are written as "shapefactor": the constrained parameters name their constraint distribution in + // "constraints", the others are free. + for (RooAbsReal *g : paramHistFuncParamsInBinOrder(*phf)) { sys.parameters.push_back(static_cast(g)); - RooAbsPdf *constraint = nullptr; - if (isshapesys) { - constraint = findConstraint(g); - if (!constraint) - constraint = ws->pdf(constraintName(g->GetName())); - if (!constraint && !g->isConstant()) { - RooJSONFactoryWSTool::error("cannot find constraint for " + std::string(g->GetName())); - } - } else if (!isshapefactor) { - RooJSONFactoryWSTool::error("unknown type of shapesys " + std::string(phf->GetName())); - } - if (!constraint) { - sys.constraints.push_back(0.0); - sys.constraintPdfs.push_back(nullptr); - } else { - sys.constraints.push_back(constraintRelError(*constraint, *g)); - sys.constraintPdfs.push_back(constraint); - } + RooAbsPdf *constraint = findConstraint(g); + if (!constraint) + constraint = ws->pdf(constraintName(g->GetName())); + sys.constraintPdfs.push_back(constraint); } sample.shapesys.emplace_back(std::move(sys)); } @@ -1687,7 +1924,7 @@ void ensureUniqueModifiers(const Channel &channel, const Sample &sample) for (const auto &modifier : sample.histosys) add("histosys", modifier.name); for (const auto &modifier : sample.shapesys) - add("shapesys", modifier.name); + add("shapefactor", modifier.name); for (const auto &modifier : sample.otherElements) add("custom", modifier.name); for (const auto &modifier : sample.tmpElements) @@ -1751,14 +1988,13 @@ std::optional defaultInterpolation(const Channel &channel) return best->first; } -bool exportChannel(RooJSONFactoryWSTool *tool, const Channel &channel, JSONNode &elem) +bool exportChannel(RooJSONFactoryWSTool *tool, const Channel &channel, ConstraintNames const &constraintNames, + JSONNode &elem) { - // Write the constraint reference for any modifier that supports an - // external Gaussian/Poisson/etc. constraint. - auto writeConstraint = [](JSONNode &mod, auto const &sys) { - if (sys.constraint) { - mod["constraint"] << sys.constraint->GetName(); - } + // Write the parameter and the name of the constraint distribution of a modifier with a single parameter. + auto writeParameter = [&](JSONNode &mod, RooAbsReal const ¶m, RooAbsPdf const *constraint) { + mod["parameters"].set_seq().append_child() << param.GetName(); + exportConstraintNames(mod, {constraint}, constraintNames); }; auto addModifier = [](JSONNode &modifiers, std::string const &name, const char *type) -> JSONNode & { auto &mod = modifiers.append_child(); @@ -1769,10 +2005,6 @@ bool exportChannel(RooJSONFactoryWSTool *tool, const Channel &channel, JSONNode }; elem["type"] << "histfactory_dist"; - // The constraint type of the staterror is also written on the channel level, where older RooFit versions read it. - if (!channel.statErrorConstraintType.empty()) { - elem[::Literals::staterror].set_map()["constraint_type"] << channel.statErrorConstraintType; - } const auto channelDefaultInterpolation = defaultInterpolation(channel); if (channelDefaultInterpolation) { writeInterpolation(elem["default_interpolation"], *channelDefaultInterpolation); @@ -1790,21 +2022,16 @@ bool exportChannel(RooJSONFactoryWSTool *tool, const Channel &channel, JSONNode auto &mod = modifiers.append_child(); mod.set_map(); mod["name"] << nf.name; - mod["parameter"] << nf.param->GetName(); + writeParameter(mod, *nf.param, nf.constraint); mod["type"] << "normfactor"; - if (nf.constraint) { - mod["constraint"] << nf.constraint->GetName(); - tool->queueExport(*nf.constraint); - } } for (const auto &sys : sample.normsys) { auto &mod = addModifier(modifiers, sys.name, "normsys"); - mod["parameter"] << sys.param->GetName(); + writeParameter(mod, *sys.param, sys.constraint); if (!channelDefaultInterpolation || sys.interpolation != *channelDefaultInterpolation) { writeInterpolation(mod["interpolation"], sys.interpolation); } - writeConstraint(mod, sys); auto &data = mod["data"].set_map(); data["lo"] << sys.low; data["hi"] << sys.high; @@ -1812,11 +2039,10 @@ bool exportChannel(RooJSONFactoryWSTool *tool, const Channel &channel, JSONNode for (const auto &sys : sample.histosys) { auto &mod = addModifier(modifiers, sys.name, "histosys"); - mod["parameter"] << sys.param->GetName(); + writeParameter(mod, *sys.param, sys.constraint); if (!channelDefaultInterpolation || sys.interpolation != *channelDefaultInterpolation) { writeInterpolation(mod["interpolation"], sys.interpolation); } - writeConstraint(mod, sys); auto &data = mod["data"].set_map(); if (channel.nBins != sys.low.size() || channel.nBins != sys.high.size()) { RooJSONFactoryWSTool::error("inconsistent binning: " + std::to_string(channel.nBins) + @@ -1828,20 +2054,10 @@ bool exportChannel(RooJSONFactoryWSTool *tool, const Channel &channel, JSONNode } for (const auto &sys : sample.shapesys) { - auto &mod = addModifier(modifiers, sys.name, "shapesys"); + auto &mod = addModifier(modifiers, sys.name, "shapefactor"); exportGammaParameters(mod, sys.parameters); - if (std::any_of(sys.constraintPdfs.begin(), sys.constraintPdfs.end(), - [](auto *pdf) { return pdf != nullptr; })) { - auto &constraintNames = mod["constraints"].set_seq(); - for (auto *constraint : sys.constraintPdfs) { - if (constraint) { - constraintNames.append_child() << constraint->GetName(); - } else { - constraintNames.append_child().set_null(); - } - } - } - mod["data"].set_map()["vals"].fill_seq(sys.constraints); + // The relative uncertainties are encoded in the constraint distributions. + exportConstraintNames(mod, sys.constraintPdfs, constraintNames); } for (const auto &other : sample.otherElements) { @@ -1854,9 +2070,7 @@ bool exportChannel(RooJSONFactoryWSTool *tool, const Channel &channel, JSONNode if (sample.useBarlowBeestonLight) { auto &mod = addModifier(modifiers, ::Literals::staterror, ::Literals::staterror); exportGammaParameters(mod, sample.staterrorParameters); - if (!channel.statErrorConstraintType.empty()) { - mod["constraint"] << channel.statErrorConstraintType; - } + exportConstraintNames(mod, sample.staterrorConstraints, constraintNames); } if (!observablesWritten) { @@ -1987,36 +2201,59 @@ bool tryExportHistFactory(RooJSONFactoryWSTool *tool, const std::string &pdfname // stat error handling configureStatError(channel); - auto lostConstraints = findLostConstraints(channel, constraints); - // Export all the lost constraints - for (const auto *constraint : lostConstraints) { - RooJSONFactoryWSTool::warning( - "losing constraint term '" + std::string(constraint->GetName()) + - "', implicit constraints are not supported by HS3 yet! The term will appear in the HS3 file, but will not be " - "picked up when creating a likelihood from it! You will have to add it manually as an external constraint."); - tool->queueExport(*constraint); - } - - // Export all the regular modifiers - auto queueConstraints = [&](auto const &modifiers) { - for (auto &modifier : modifiers) { - if (modifier.constraint) { - tool->queueExport(*modifier.constraint); + // Export the constraints of the modifiers, in the order in which the modifiers are written. The terms under + // which they are written are recorded for the likelihoods of the models that contain this channel. + ConstraintNames constraintNames; + std::vector terms; + auto exportConstraintOf = [&](RooAbsReal const *param, RooAbsPdf const *constraint) { + if (constraint && constraintNames.find(constraint) == constraintNames.end()) { + // HistFactory builds the Gaussian constraints with the parameter as x and the global observable as mean. + // The observable is the side that does not depend on the parameter; recording it makes the Gaussian be + // written with the observable as x. + if (auto const *gauss = dynamic_cast(constraint)) { + RooAbsReal const &observable = gauss->getX().dependsOn(*param) ? gauss->getMean() : gauss->getX(); + if (!isLiteralConstVar(observable)) { + tool->setConstraintObservable(*constraint, observable.GetName()); + } + } + const std::string term = exportConstraint(*tool, *constraint); + constraintNames[constraint] = term; + if (std::find(terms.begin(), terms.end(), term) == terms.end()) { + terms.push_back(term); } } }; for (const auto &sample : channel.samples) { - queueConstraints(sample.normfactors); - queueConstraints(sample.normsys); - queueConstraints(sample.histosys); + for (auto &modifier : sample.normfactors) { + exportConstraintOf(modifier.param, modifier.constraint); + } + for (auto &modifier : sample.normsys) { + exportConstraintOf(modifier.param, modifier.constraint); + } + for (auto &modifier : sample.histosys) { + exportConstraintOf(modifier.param, modifier.constraint); + } for (auto &modifier : sample.shapesys) { - for (auto *constraint : modifier.constraintPdfs) { - if (constraint) { - tool->queueExport(*constraint); - } + for (std::size_t i = 0; i < modifier.constraintPdfs.size(); ++i) { + exportConstraintOf(modifier.parameters[i], modifier.constraintPdfs[i]); } } + for (std::size_t i = 0; i < sample.staterrorConstraints.size(); ++i) { + exportConstraintOf(sample.staterrorParameters[i], sample.staterrorConstraints[i]); + } + } + + // Constraint terms that are not attached to any modifier parameter are written on their own; the likelihood + // lists them with the other constraint terms. + for (const auto *constraint : findLostConstraints(channel, constraints)) { + oocoutI(nullptr, IO) << "RooFitHS3: the constraint term '" << constraint->GetName() << "' of '" << pdfname + << "' is not attached to any modifier parameter; it is written on its own and listed in " + "the likelihood" + << std::endl; + tool->queueExport(*constraint); + terms.push_back(constraint->GetName()); } + tool->recordConstraintTerms(pdfname, terms); // Export all the custom modifiers for (const auto &sample : channel.samples) { @@ -2038,7 +2275,7 @@ bool tryExportHistFactory(RooJSONFactoryWSTool *tool, const std::string &pdfname } } - return exportChannel(tool, channel, elem); + return exportChannel(tool, channel, constraintNames, elem); } class HistFactoryStreamer_ProdPdf : public RooFit::JSONIO::Exporter { diff --git a/roofit/hs3/src/JSONFactories_RooFitCore.cxx b/roofit/hs3/src/JSONFactories_RooFitCore.cxx index ce7102885bc66..b6473b334e87a 100644 --- a/roofit/hs3/src/JSONFactories_RooFitCore.cxx +++ b/roofit/hs3/src/JSONFactories_RooFitCore.cxx @@ -17,6 +17,8 @@ #include #include #include +#include +#include #include #include #include @@ -26,6 +28,7 @@ #include #include #include +#include #include #include #include @@ -57,7 +60,9 @@ #include "JSONIOUtils.h" +#include #include +#include #include "static_execute.h" @@ -66,7 +71,9 @@ #include #include #include +#include #include +#include #include #include @@ -402,6 +409,203 @@ bool importProdPdf(RooJSONFactoryWSTool *tool, const JSONNode &p) return true; } +// The name of the pdf that represents the component of an array distribution over the observable `observable`. +std::string arrayComponentName(std::string const &observable) +{ + return observable + "_constraint"; +} + +// Copy the value of a JSON node, keeping numbers and booleans as such. +void copyValue(JSONNode &target, JSONNode const &source) +{ + if (source.is_number()) { + target << source.val_double(); + return; + } + const std::string value = source.val(); + if (value == "true" || value == "false") { + target << (value == "true"); + } else { + target << value; + } +} + +// An HS3 `array_dist` is a product of one-dimensional distributions of one type (`component_type`), each over its own +// observable `x`. RooFit has no dedicated class for it: every component is imported as an individual pdf of the +// component type, named after its observable, and a RooProdPdf with the attribute `arrayDistTag` that carries the +// name of the array holds them together, so that the array can be written back as `array_dist`. HistFactory modifiers +// that name an array as the constraint of a parameter pick their component out of it (see +// JSONFactories_HistFactory.cxx). +bool importArrayDist(RooJSONFactoryWSTool *tool, const JSONNode &p) +{ + const std::string name = RooJSONFactoryWSTool::name(p); + auto error = [&name](std::string const &message) { + RooJSONFactoryWSTool::error("array_dist '" + name + "': " + message); + }; + + static const std::set supportedComponentTypes{"poisson_dist", "gaussian_dist", "normal_dist", + "lognormal_dist"}; + if (!p.has_child("component_type")) { + error("missing \"component_type\""); + } + const std::string componentType = p["component_type"].val(); + if (supportedComponentTypes.find(componentType) == supportedComponentTypes.end()) { + error("the component_type '" + componentType + + "' is not supported, supported are poisson_dist, gaussian_dist and lognormal_dist"); + } + + const JSONNode *defaults = p.find("default_parameters"); + if (defaults) { + if (!defaults->is_map()) { + error("\"default_parameters\" must be a struct"); + } + if (defaults->has_child("x")) { + error("\"x\" must be given per component, not in \"default_parameters\""); + } + } + const JSONNode *components = p.find("components"); + if (!components || !components->is_seq() || components->num_children() == 0) { + error("\"components\" must be a non-empty list"); + } + + RooWorkspace &ws = *tool->workspace(); + + // The components are imported through the importer of the component type, from a node per component that + // combines its own fields with the defaults. + auto tree = RooFit::Detail::JSONTree::create(); + JSONNode &nodes = tree->rootnode().set_seq(); + std::vector observables; + RooArgList pdfs; + for (JSONNode const &component : components->children()) { + const std::string context = "component " + std::to_string(observables.size()); + if (!component.is_map()) { + error(context + " must be a struct"); + } + const JSONNode *xNode = component.find("x"); + if (!xNode) { + error(context + " has no \"x\""); + } + if (xNode->is_number() || !xNode->has_val()) { + error(context + ": \"x\" must be the name of the observable"); + } + const std::string x = xNode->val(); + if (std::find(observables.begin(), observables.end(), x) != observables.end()) { + error("the observable '" + x + "' is used by two components"); + } + observables.push_back(x); + + const std::string pdfName = arrayComponentName(x); + if (ws.arg(pdfName)) { + error("cannot create the component for '" + x + "' as '" + pdfName + + "', the workspace already has an object of that name"); + } + JSONNode &node = nodes.append_child().set_map(); + node["name"] << pdfName; + node["type"] << componentType; + if (defaults) { + for (JSONNode const &field : defaults->children()) { + copyValue(node[field.key()], field); + } + } + for (JSONNode const &field : component.children()) { + copyValue(node[field.key()], field); + } + // RooFit's "integer" flag of poisson_dist (not part of HS3) defaults to the continued Poisson distribution, + // which is what HistFactory constraints are and which gives the usual value at integer x. + if (componentType == "poisson_dist" && !node.has_child("integer")) { + node["integer"] << false; + } + tool->importFunction(node, true); + RooAbsPdf *pdf = ws.pdf(pdfName); + if (!pdf) { + error(context + " did not produce a distribution"); + } + pdfs.add(*pdf); + } + + // The parameters of a component must not depend on any observable of the array, otherwise the product would not + // be normalized. + for (std::size_t i = 0; i < observables.size(); ++i) { + RooAbsArg const *x = ws.arg(observables[i]); + if (!x) { + error("the observable '" + observables[i] + "' is not a parameter or function"); + } + for (std::size_t k = 0; k < observables.size(); ++k) { + if (k != i && pdfs[i].dependsOn(*ws.arg(observables[k]))) { + error("component " + std::to_string(i) + " depends on the observable '" + observables[k] + + "' of component " + std::to_string(k)); + } + } + for (RooAbsArg const *server : pdfs[i].servers()) { + if (server != x && server->dependsOn(*x)) { + error("the parameters of component " + std::to_string(i) + " depend on its own observable '" + + observables[i] + "'"); + } + } + } + + RooProdPdf array{name.c_str(), name.c_str(), pdfs}; + array.setAttribute(arrayDistTag); + tool->wsImport(array); + return true; +} + +// One field of a component of an array distribution: a name, a number or a flag. +using ArrayFieldValue = std::variant; + +struct ArrayComponent { + std::string type; + std::string x; + std::vector> fields; +}; + +ArrayFieldValue nameOrNumber(RooAbsReal const &arg) +{ + if (isLiteralConstVar(arg)) { + return static_cast(arg).getVal(); + } + return std::string{arg.GetName()}; +} + +// The type, observable and fields of `pdf` as a component of an array distribution, with the fields of the standalone +// exporters, or std::nullopt if `pdf` cannot be a component. A Gaussian whose mean is the recorded `observable` is +// written with x and mean swapped. +std::optional arrayComponentOf(RooAbsPdf const &pdf, std::string const *observable = nullptr) +{ + ArrayComponent out; + RooAbsReal const *x = nullptr; + if (auto const *poisson = dynamic_cast(&pdf)) { + out.type = "poisson_dist"; + x = &poisson->getX(); + out.fields = {{"mean", nameOrNumber(poisson->getMean())}, {"integer", !poisson->getNoRounding()}}; + } else if (auto const *gauss = dynamic_cast(&pdf)) { + out.type = "gaussian_dist"; + x = &gauss->getX(); + RooAbsReal const *mean = &gauss->getMean(); + // A Gaussian constraint built by HistFactory has the parameter as x and the global observable as mean. The + // Gaussian is symmetric in the two, so it is written with the observable as x where the observable is known. + if (observable && *observable == mean->GetName() && !isLiteralConstVar(*mean)) { + std::swap(x, mean); + } + out.fields = {{"mean", nameOrNumber(*mean)}, {"sigma", nameOrNumber(gauss->getSigma())}}; + } else if (auto const *lognormal = dynamic_cast(&pdf)) { + if (!lognormal->useStandardParametrization()) { + return std::nullopt; + } + out.type = "lognormal_dist"; + x = &lognormal->getX(); + out.fields = {{"mu", nameOrNumber(lognormal->getMedian())}, {"sigma", nameOrNumber(lognormal->getShapeK())}}; + } else { + return std::nullopt; + } + // The observable of a component must be named: the observables are pairwise distinct and are paired with data. + if (isLiteralConstVar(*x)) { + return std::nullopt; + } + out.x = x->GetName(); + return out; +} + bool importAddition(RooJSONFactoryWSTool *tool, const JSONNode &p) { std::string name(RooJSONFactoryWSTool::name(p)); @@ -789,10 +993,16 @@ bool importParamHistFunc(RooJSONFactoryWSTool *tool, const JSONNode &p) << ". Using default binning (uniform; nbins=100). If needed, export the Workspace to JSON with a newer " << "Root version that supports custom ParamHistFunc binnings(>=6.38.00)." << std::endl; RooJSONFactoryWSTool::warning(ss.str()); - tool->wsEmplace(name, varList, tool->requestArgList(p, "parameters")); + checkParamHistFuncOrderIsDefined(tool->hs3VersionInput(), varList, "the ParamHistFunc '" + name + "'"); + tool->wsEmplace(name, varList, + toParamHistFuncOrder(tool->requestArgList(p, "parameters"), varList)); return true; } - tool->wsEmplace(name, readBinning(p, varList), tool->requestArgList(p, "parameters")); + RooArgList binning = readBinning(p, varList); + checkParamHistFuncOrderIsDefined(tool->hs3VersionInput(), binning, "the ParamHistFunc '" + name + "'"); + // The parameters are given in the order of HS3 (last axis fastest); ParamHistFunc expects TH1 order. + tool->wsEmplace(name, binning, + toParamHistFuncOrder(tool->requestArgList(p, "parameters"), binning)); return true; } @@ -869,9 +1079,152 @@ class DeprecatedKeyImporter : public RooFit::JSONIO::Importer { mutable bool _warned = false; }; +namespace CrystalBall { + +constexpr const char *tailKeys[] = {"alpha_L", "n_L", "alpha_R", "n_R"}; + +bool hasAny(const JSONNode &p, std::initializer_list keys) +{ + for (const char *key : keys) { + if (p.has_child(key)) { + return true; + } + } + return false; +} + +bool hasDoubleSidedKeys(const JSONNode &p) +{ + return hasAny(p, {"alpha_L", "n_L", "alpha_R", "n_R", "sigma_L", "sigma_R"}); +} + +// Create a RooCrystalBall with two tails. With `symmetricShorthand`, the tails may also be given as "alpha" and "n" +// for both sides, as in the "crystalball_dist" of HS3 versions before 0.3. +void importDoubleSided(RooJSONFactoryWSTool *tool, const JSONNode &p, std::string const &type, bool symmetricShorthand) +{ + const std::string name(RooJSONFactoryWSTool::name(p)); + const std::string context = type + " '" + name + "'"; + + const bool asymmetricTails = hasAny(p, {"alpha_L", "n_L", "alpha_R", "n_R"}); + const bool symmetricTails = hasAny(p, {"alpha", "n"}); + if (symmetricTails && !symmetricShorthand) { + RooJSONFactoryWSTool::error(context + " must not have \"alpha\" or \"n\", the tails are given by \"alpha_L\", " + + "\"n_L\", \"alpha_R\" and \"n_R\""); + } + if (symmetricTails && asymmetricTails) { + RooJSONFactoryWSTool::error(context + + " must not combine \"alpha\"/\"n\" with \"alpha_L\"/\"n_L\"/\"alpha_R\"/\"n_R\""); + } + if (asymmetricTails || !symmetricTails) { + for (const char *key : tailKeys) { + if (!p.has_child(key)) { + RooJSONFactoryWSTool::error(context + " is missing \"" + key + + "\": \"alpha_L\", \"n_L\", \"alpha_R\" and \"n_R\" are all required"); + } + } + } else if (!p.has_child("alpha") || !p.has_child("n")) { + RooJSONFactoryWSTool::error(context + " needs both \"alpha\" and \"n\""); + } + + const bool asymmetricWidths = hasAny(p, {"sigma_L", "sigma_R"}); + if (asymmetricWidths && p.has_child("sigma")) { + RooJSONFactoryWSTool::error(context + " must not combine \"sigma\" with \"sigma_L\"/\"sigma_R\""); + } + if (asymmetricWidths && !(p.has_child("sigma_L") && p.has_child("sigma_R"))) { + RooJSONFactoryWSTool::error(context + " needs both \"sigma_L\" and \"sigma_R\""); + } + if (!asymmetricWidths && !p.has_child("sigma")) { + RooJSONFactoryWSTool::error(context + " needs either \"sigma\" or both \"sigma_L\" and \"sigma_R\""); + } + + RooAbsReal *m = tool->requestArg(p, "m"); + RooAbsReal *m0 = tool->requestArg(p, "m0"); + RooAbsReal *sigmaL = tool->requestArg(p, asymmetricWidths ? "sigma_L" : "sigma"); + RooAbsReal *sigmaR = tool->requestArg(p, asymmetricWidths ? "sigma_R" : "sigma"); + RooAbsReal *alphaL = tool->requestArg(p, symmetricTails ? "alpha" : "alpha_L"); + RooAbsReal *nL = tool->requestArg(p, symmetricTails ? "n" : "n_L"); + RooAbsReal *alphaR = tool->requestArg(p, symmetricTails ? "alpha" : "alpha_R"); + RooAbsReal *nR = tool->requestArg(p, symmetricTails ? "n" : "n_R"); + tool->wsEmplace(name, *m, *m0, *sigmaL, *sigmaR, *alphaL, *nL, *alphaR, *nR); +} + +} // namespace CrystalBall + +// "crystalball_dist" is the single-sided Crystal Ball, which is left to the factory expression (RooCBShape). Before +// HS3 0.3, it was also the double-sided Crystal Ball, recognizable by the "_L"/"_R" fields. Such documents are still +// read, with a deprecation warning. +class CrystalBallImporter : public RooFit::JSONIO::Importer { +public: + bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override + { + if (!CrystalBall::hasDoubleSidedKeys(p)) { + return false; + } + if (!_warned) { + RooJSONFactoryWSTool::warning("RooFitHS3: \"crystalball_dist\" with the fields \"alpha_L\", \"n_L\", " + "\"alpha_R\", \"n_R\", \"sigma_L\" or \"sigma_R\" is deprecated, please use " + "\"crystalball_doublesided_dist\" for the double-sided Crystal Ball."); + _warned = true; + } + CrystalBall::importDoubleSided(tool, p, "crystalball_dist", /*symmetricShorthand=*/true); + return true; + } + +private: + mutable bool _warned = false; +}; + +bool importCrystalBallDoubleSided(RooJSONFactoryWSTool *tool, const JSONNode &p) +{ + CrystalBall::importDoubleSided(tool, p, "crystalball_doublesided_dist", /*symmetricShorthand=*/false); + return true; +} + /////////////////////////////////////////////////////////////////////////////////////////////////////// // specialized exporter implementations /////////////////////////////////////////////////////////////////////////////////////////////////////// +// Write a reference to `arg`: the value if it is a literal constant, and the name otherwise. +void writeReference(JSONNode &node, RooAbsArg const &arg) +{ + auto const *constVar = dynamic_cast(&arg); + std::istringstream stream(arg.GetName()); + double value = 0.0; + if (constVar && (stream >> value) && stream.eof()) { + node << constVar->getVal(); + } else { + node << arg.GetName(); + } +} + +// A RooCrystalBall with one tail is the single-sided "crystalball_dist" (read back as RooCBShape, which follows the +// same convention for the sign of alpha). With two tails, it is a "crystalball_doublesided_dist", whose tails are +// always written as "alpha_L"/"n_L"/"alpha_R"/"n_R": there is no shorthand for symmetric tails. +bool exportCrystalBall(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &elem, std::string const &key) +{ + auto const *pdf = static_cast(func); + writeReference(elem["m"], pdf->x()); + writeReference(elem["m0"], pdf->x0()); + if (!pdf->alphaR() || !pdf->nR()) { + elem["type"] << "crystalball_dist"; + writeReference(elem["sigma"], pdf->sigmaL()); + writeReference(elem["alpha"], pdf->alphaL()); + writeReference(elem["n"], pdf->nL()); + return true; + } + elem["type"] << key; + if (&pdf->sigmaL() == &pdf->sigmaR()) { + writeReference(elem["sigma"], pdf->sigmaL()); + } else { + writeReference(elem["sigma_L"], pdf->sigmaL()); + writeReference(elem["sigma_R"], pdf->sigmaR()); + } + writeReference(elem["alpha_L"], pdf->alphaL()); + writeReference(elem["n_L"], pdf->nL()); + writeReference(elem["alpha_R"], *pdf->alphaR()); + writeReference(elem["n_R"], *pdf->nR()); + return true; +} + template bool exportAddPdf(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &elem, std::string const &key) { @@ -1205,6 +1558,25 @@ bool exportGaussModel(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &e return true; } +// A Gaussian constraint whose recorded observable (see RooJSONFactoryWSTool::setConstraintObservable()) is its +// mean is written with x and mean swapped, so that the observable is x. Other Gaussians use the generic export. +bool exportGaussianConstraint(RooJSONFactoryWSTool *tool, const RooAbsArg *func, JSONNode &elem, std::string const &key) +{ + auto const *gauss = static_cast(func); + std::string const *observable = tool->constraintObservable(gauss->GetName()); + if (!observable || *observable != gauss->getMean().GetName() || isLiteralConstVar(gauss->getMean())) { + return false; + } + auto write = [](JSONNode &node, ArrayFieldValue const &value) { + std::visit([&node](auto const &v) { node << v; }, value); + }; + elem["type"] << key; + elem["x"] << gauss->getMean().GetName(); + write(elem["mean"], nameOrNumber(gauss->getX())); + write(elem["sigma"], nameOrNumber(gauss->getSigma())); + return true; +} + bool exportLogNormal(RooJSONFactoryWSTool *tool, const RooAbsArg *func, JSONNode &elem, std::string const &key) { auto *pdf = static_cast(func); @@ -1329,7 +1701,11 @@ bool exportParamHistFunc(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode auto *pdf = static_cast(func); elem["type"] << key; RooJSONFactoryWSTool::fillSeq(elem["variables"], pdf->dataVars()); - RooJSONFactoryWSTool::fillSeq(elem["parameters"], pdf->paramList()); + // One parameter per bin, in the order of HS3 (last axis fastest), not in ParamHistFunc's internal TH1 order. + auto ¶msNode = elem["parameters"].set_seq(); + for (RooAbsReal *param : paramHistFuncParamsInBinOrder(*pdf)) { + paramsNode.append_child() << param->GetName(); + } auto &observablesNode = elem["axes"].set_seq(); // axes have to be ordered to get consistent bin indices for (auto *var : static_range_cast(pdf->dataVars())) { @@ -1459,11 +1835,14 @@ void registerExporter(TClass const *cl, std::string key, bool topPriority = true STATIC_EXECUTE([]() { registerImporter("density_function_dist"); registerImporter("rate_extended_dist"); + RooFit::JSONIO::registerImporter("crystalball_dist", std::make_unique()); + registerImporter("crystalball_doublesided_dist"); RooFit::JSONIO::registerImporter("ARGUS_dist", std::make_unique("ARGUS_dist", "argus_dist")); RooFit::JSONIO::registerImporter("binsampling", std::make_unique("binsampling", "binsampling_dist")); registerImporter("product", false); registerImporter("product_dist", false); + registerImporter("array_dist", false); registerImporter("sum", false); registerImporter("mixture_dist", false); registerImporter("mixture_resolution_model", false); @@ -1500,6 +1879,7 @@ STATIC_EXECUTE([]() { registerExporter>(RooAddPdf::Class(), "mixture_dist", false); registerExporter>(RooAddModel::Class(), "mixture_resolution_model", false); registerExporter(RooBinSamplingPdf::Class(), "binsampling_dist", false); + registerExporter(RooCrystalBall::Class(), "crystalball_doublesided_dist", false); registerExporter(RooBinWidthFunction::Class(), "binvolume", false); registerExporter>(RooLegacyExpPoly::Class(), "legacy_exp_poly_dist", false); registerExporter(RooExponential::Class(), "exponential_dist", false); @@ -1508,6 +1888,7 @@ STATIC_EXECUTE([]() { registerExporter>(RooHistFunc::Class(), "step", false); registerExporter>(RooHistPdf::Class(), "histogram_dist", false); registerExporter(RooLognormal::Class(), "lognormal_dist", false); + registerExporter(RooGaussian::Class(), "gaussian_dist", true); registerExporter(RooMultiVarGaussian::Class(), "multivariate_normal_dist", false); registerExporter(RooPoisson::Class(), "poisson_dist", false); registerExporter(RooDecay::Class(), "decay_dist", false, /*autoExportDependants=*/false); @@ -1527,3 +1908,80 @@ STATIC_EXECUTE([]() { }); } // namespace + +std::optional arrayComponentInfo(RooAbsPdf const &pdf, std::string const *observable) +{ + auto component = arrayComponentOf(pdf, observable); + if (!component) { + return std::nullopt; + } + return ArrayComponentInfo{component->type, component->x}; +} + +void writeArrayDist(JSONNode &elem, std::string const &name, std::vector const &components, + std::map const &observableOf) +{ + auto error = [&name](std::string const &message) { + RooJSONFactoryWSTool::error("cannot write the array distribution '" + name + "': " + message); + }; + if (components.empty()) { + error("it has no components"); + } + std::vector comps; + std::set observables; + for (RooAbsPdf const *pdf : components) { + auto recorded = observableOf.find(pdf->GetName()); + auto comp = arrayComponentOf(*pdf, recorded != observableOf.end() ? &recorded->second : nullptr); + if (!comp) { + error("'" + std::string(pdf->GetName()) + + "' cannot be a component, only RooPoisson, RooGaussian and RooLognormal (standard parametrization) " + "over a named observable are supported"); + } + if (!comps.empty() && comp->type != comps.front().type) { + error("the components '" + std::string(components.front()->GetName()) + "' (" + comps.front().type + + ") and '" + pdf->GetName() + "' (" + comp->type + ") have different types"); + } + if (!observables.insert(comp->x).second) { + error("two components have the observable '" + comp->x + "'"); + } + comps.push_back(std::move(*comp)); + } + + elem["type"] << "array_dist"; + elem["component_type"] << comps.front().type; + + auto write = [](JSONNode &node, ArrayFieldValue const &value) { + std::visit([&node](auto const &v) { node << v; }, value); + }; + + // A field with the same value in every component is written once in "default_parameters". + std::set defaultKeys; + if (comps.size() > 1) { + for (std::size_t i = 0; i < comps.front().fields.size(); ++i) { + auto const &field = comps.front().fields[i]; + if (std::all_of(comps.begin(), comps.end(), + [&](ArrayComponent const &comp) { return comp.fields[i].second == field.second; })) { + defaultKeys.insert(field.first); + } + } + } + if (!defaultKeys.empty()) { + JSONNode &defaults = elem["default_parameters"].set_map(); + for (auto const &field : comps.front().fields) { + if (defaultKeys.find(field.first) != defaultKeys.end()) { + write(defaults[field.first], field.second); + } + } + } + + JSONNode &seq = elem["components"].set_seq(); + for (ArrayComponent const &comp : comps) { + JSONNode &node = seq.append_child().set_map(); + node["x"] << comp.x; + for (auto const &field : comp.fields) { + if (defaultKeys.find(field.first) == defaultKeys.end()) { + write(node[field.first], field.second); + } + } + } +} diff --git a/roofit/hs3/src/JSONFactories_RooFitCore.cxx~ b/roofit/hs3/src/JSONFactories_RooFitCore.cxx~ new file mode 100644 index 0000000000000..df4d823d6260c --- /dev/null +++ b/roofit/hs3/src/JSONFactories_RooFitCore.cxx~ @@ -0,0 +1,1350 @@ +/* + * Project: RooFit + * Authors: + * Carsten D. Burgard, DESY/ATLAS, Dec 2021 + * + * Copyright (c) 2022, CERN + * + * Redistribution and use in source and binary forms, + * with or without modification, are permitted according to the terms + * listed in LICENSE (http://roofit.sourceforge.net/license.txt) + */ + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include "JSONIOUtils.h" + +#include "static_execute.h" + +#include +#include + +using RooFit::Detail::JSONNode; + +/////////////////////////////////////////////////////////////////////////////////////////////////////// +// individually implemented importers +/////////////////////////////////////////////////////////////////////////////////////////////////////// + +namespace { +/** + * Extracts arguments from a mathematical expression. + * + * This function takes a string representing a mathematical + * expression and extracts the arguments from it. The arguments are + * defined as sequences of characters that do not contain digits, + * spaces, or parentheses, and that start with a letter. Function + * calls such as "exp( ... )", identified as being followed by an + * opening parenthesis, are not treated as arguments. The extracted + * arguments are returned as a vector of strings. + * + * @param expr A string representing a mathematical expression. + * @return A set of unique strings representing the extracted arguments. + */ +std::set extractArguments(std::string expr) +{ + // Get rid of whitespaces + expr.erase(std::remove_if(expr.begin(), expr.end(), [](unsigned char c) { return std::isspace(c); }), expr.end()); + + std::set arguments; + size_t startidx = expr.size(); + for (size_t i = 0; i < expr.size(); ++i) { + if (startidx >= expr.size()) { + if (isalpha(expr[i])) { + startidx = i; + // check this character is not part of scientific notation, e.g. 2e-5 + if (TFormula::IsScientificNotation(expr, i)) { + // if it is, we ignore this character + startidx = expr.size(); + } + } + } else { + if (!isdigit(expr[i]) && !isalpha(expr[i]) && expr[i] != '_') { + if (expr[i] == '(') { + startidx = expr.size(); + continue; + } + std::string arg(expr.substr(startidx, i - startidx)); + startidx = expr.size(); + arguments.insert(arg); + } + } + } + if (startidx < expr.size()) { + arguments.insert(expr.substr(startidx)); + } + return arguments; +} + +template +class RooFormulaArgFactory : public RooFit::JSONIO::Importer { +public: + bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override + { + std::string name(RooJSONFactoryWSTool::name(p)); + if (!p.has_child("expression")) { + RooJSONFactoryWSTool::error("no expression given for '" + name + "'"); + } + TString formula(p["expression"].val()); + RooArgList dependents; + for (const auto &d : extractArguments(formula.Data())) { + dependents.add(*tool->request(d, name)); + } + tool->wsImport(RooArg_t{name.c_str(), formula, dependents}); + return true; + } +}; + +// Fast-path importers for RooProduct, RooAddition, and RooProdPdf that +// bypass the generic factory-expression mechanism. The default path +// generates a string expression and passes it to gROOT->ProcessLineFast(), +// which invokes the Cling JIT for every single call. For workspaces with +// thousands of product/sum nodes (a common shape for HistFactory models) +// that JIT cost dominates JSON import time. Constructing the RooFit object +// directly here keeps the work O(N) of cheap C++ calls. +class RooProductFactory : public RooFit::JSONIO::Importer { +public: + bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override + { + std::string name(RooJSONFactoryWSTool::name(p)); + tool->wsEmplace(name, tool->requestArgList(p, "factors")); + return true; + } +}; + +class RooProdPdfFactory : public RooFit::JSONIO::Importer { +public: + bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override + { + std::string name(RooJSONFactoryWSTool::name(p)); + tool->wsEmplace(name, tool->requestArgList(p, "factors")); + return true; + } +}; + +class RooAdditionFactory : public RooFit::JSONIO::Importer { +public: + bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override + { + std::string name(RooJSONFactoryWSTool::name(p)); + tool->wsEmplace(name, tool->requestArgList(p, "summands")); + return true; + } +}; + +class RooAddPdfFactory : public RooFit::JSONIO::Importer { +public: + bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override + { + std::string name(RooJSONFactoryWSTool::name(p)); + if (!tool->requestArgList(p, "coefficients").empty()) { + tool->wsEmplace(name, tool->requestArgList(p, "summands"), + tool->requestArgList(p, "coefficients")); + return true; + } + tool->wsEmplace(name, tool->requestArgList(p, "summands")); + return true; + } +}; + +class RooAddModelFactory : public RooFit::JSONIO::Importer { +public: + bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override + { + std::string name(RooJSONFactoryWSTool::name(p)); + tool->wsEmplace(name, tool->requestArgList(p, "summands"), + tool->requestArgList(p, "coefficients")); + return true; + } +}; + +class RooBinWidthFunctionFactory : public RooFit::JSONIO::Importer { +public: + bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override + { + std::string name(RooJSONFactoryWSTool::name(p)); + RooHistFunc *hf = static_cast(tool->request(p["histogram"].val(), name)); + tool->wsEmplace(name, *hf, p["divideByBinWidth"].val_bool()); + return true; + } +}; + +class RooBinSamplingPdfFactory : public RooFit::JSONIO::Importer { +public: + bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override + { + std::string name(RooJSONFactoryWSTool::name(p)); + + RooAbsPdf *pdf = tool->requestArg(p, "pdf"); + RooRealVar *obs = tool->requestArg(p, "observable"); + + if (!pdf->dependsOn(*obs)) { + RooJSONFactoryWSTool::error(std::string("pdf '") + pdf->GetName() + "' does not depend on observable '" + + obs->GetName() + "' as indicated by parent RooBinSamplingPdf '" + name + + "', please check!"); + } + + if (!p.has_child("epsilon")) { + RooJSONFactoryWSTool::error("no epsilon given in '" + name + "'"); + } + double epsilon(p["epsilon"].val_double()); + + tool->wsEmplace(name, *obs, *pdf, epsilon); + + return true; + } +}; + +class RooRealSumPdfFactory : public RooFit::JSONIO::Importer { +public: + bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override + { + std::string name(RooJSONFactoryWSTool::name(p)); + + bool extended = false; + if (p.has_child("extended") && p["extended"].val_bool()) { + extended = true; + } + tool->wsEmplace(name, tool->requestArgList(p, "samples"), + tool->requestArgList(p, "coefficients"), extended); + return true; + } +}; + +class RooRealSumFuncFactory : public RooFit::JSONIO::Importer { +public: + bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override + { + std::string name(RooJSONFactoryWSTool::name(p)); + tool->wsEmplace(name, tool->requestArgList(p, "samples"), + tool->requestArgList(p, "coefficients")); + return true; + } +}; +template +class RooPolynomialFactory : 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 RooPoissonFactory : public RooFit::JSONIO::Importer { +public: + bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override + { + std::string name(RooJSONFactoryWSTool::name(p)); + RooAbsReal *x = tool->requestArg(p, "x"); + RooAbsReal *mean = tool->requestArg(p, "mean"); + tool->wsEmplace(name, *x, *mean, !p["integer"].val_bool()); + return true; + } +}; + +class RooDecayFactory : public RooFit::JSONIO::Importer { +public: + bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override + { + std::string name(RooJSONFactoryWSTool::name(p)); + RooRealVar *t = tool->requestArg(p, "t"); + RooAbsReal *tau = tool->requestArg(p, "tau"); + RooResolutionModel *model = dynamic_cast(tool->requestArg(p, "resolutionModel")); + RooDecay::DecayType decayType = static_cast(p["decayType"].val_int()); + tool->wsEmplace(name, *t, *tau, *model, decayType); + return true; + } +}; + +class RooTruthModelFactory : public RooFit::JSONIO::Importer { +public: + bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override + { + std::string name(RooJSONFactoryWSTool::name(p)); + RooRealVar *x = tool->requestArg(p, "x"); + tool->wsEmplace(name, *x); + return true; + } +}; + +class RooGaussModelFactory : public RooFit::JSONIO::Importer { +public: + bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override + { + std::string name(RooJSONFactoryWSTool::name(p)); + RooRealVar *x = tool->requestArg(p, "x"); + RooRealVar *mean = tool->requestArg(p, "mean"); + RooRealVar *sigma = tool->requestArg(p, "sigma"); + tool->wsEmplace(name, *x, *mean, *sigma); + return true; + } +}; + +class RooRealIntegralFactory : public RooFit::JSONIO::Importer { +public: + bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override + { + std::string name(RooJSONFactoryWSTool::name(p)); + RooAbsReal *func = tool->requestArg(p, "integrand"); + auto vars = tool->requestArgList(p, "variables"); + RooArgSet normSet; + RooArgSet const *normSetPtr = nullptr; + if (p.has_child("normalization")) { + normSet.add(tool->requestArgSet(p, "normalization")); + normSetPtr = &normSet; + } + std::string domain; + bool hasDomain = p.has_child("domain"); + if (hasDomain) { + domain = p["domain"].val(); + } + // todo: at some point, take care of integrator configurations + tool->wsEmplace(name, *func, vars, normSetPtr, static_cast(nullptr), + hasDomain ? domain.c_str() : nullptr); + return true; + } +}; + +class RooDerivativeFactory : public RooFit::JSONIO::Importer { +public: + bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override + { + std::string name(RooJSONFactoryWSTool::name(p)); + RooAbsReal *func = tool->requestArg(p, "function"); + RooRealVar *x = tool->requestArg(p, "x"); + Int_t order = p["order"].val_int(); + double eps = p["eps"].val_double(); + if (p.has_child("normalization")) { + RooArgSet normSet; + normSet.add(tool->requestArgSet(p, "normalization")); + tool->wsEmplace(name, *func, *x, normSet, order, eps); + return true; + } + tool->wsEmplace(name, *func, *x, order, eps); + return true; + } +}; + +class RooFFTConvPdfFactory : public RooFit::JSONIO::Importer { +public: + bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override + { + std::string name(RooJSONFactoryWSTool::name(p)); + RooRealVar *convVar = tool->requestArg(p, "conv_var"); + Int_t order = p["ipOrder"].val_int(); + RooAbsPdf *pdf1 = tool->requestArg(p, "pdf1"); + RooAbsPdf *pdf2 = tool->requestArg(p, "pdf2"); + if (p.has_child("conv_func")) { + RooAbsReal *convFunc = tool->requestArg(p, "conv_func"); + tool->wsEmplace(name, *convFunc, *convVar, *pdf1, *pdf2, order); + return true; + } + tool->wsEmplace(name, *convVar, *pdf1, *pdf2, order); + return true; + } +}; + +class RooExtendPdfFactory : public RooFit::JSONIO::Importer { +public: + bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override + { + std::string name(RooJSONFactoryWSTool::name(p)); + RooAbsPdf *pdf = tool->requestArg(p, "pdf"); + RooAbsReal *norm = tool->requestArg(p, "norm"); + if (p.has_child("range")) { + std::string rangeName = p["range"].val(); + tool->wsEmplace(name, *pdf, *norm, rangeName.c_str()); + return true; + } + tool->wsEmplace(name, *pdf, *norm); + return true; + } +}; + +class RooLogNormalFactory : public RooFit::JSONIO::Importer { +public: + bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override + { + std::string name(RooJSONFactoryWSTool::name(p)); + RooAbsReal *x = tool->requestArg(p, "x"); + + // Same mechanism to undo the parameter transformation as in the + // RooExponentialFactory (see comments in that class for more info). + const std::string muName = p["mu"].val(); + const std::string sigmaName = p["sigma"].val(); + const bool isTransformed = endsWith(muName, "_lognormal_log"); + const std::string suffixToRemove = isTransformed ? "_lognormal_log" : ""; + RooAbsReal *mu = tool->request(removeSuffix(muName, suffixToRemove), name); + RooAbsReal *sigma = tool->request(removeSuffix(sigmaName, suffixToRemove), name); + + tool->wsEmplace(name, *x, *mu, *sigma, !isTransformed); + + return true; + } +}; + +class RooExponentialFactory : public RooFit::JSONIO::Importer { +public: + bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override + { + std::string name(RooJSONFactoryWSTool::name(p)); + RooAbsReal *x = tool->requestArg(p, "x"); + + // If the parameter name ends with the "_exponential_inverted" suffix, + // this means that it was exported from a RooFit object where the + // parameter first needed to be transformed on export to match the HS3 + // specification. But when re-importing such a parameter, we can simply + // skip the transformation and use the original RooFit parameter without + // the suffix. + // + // A concrete example: take the following RooFit pdf in the factory language: + // + // "Exponential::exponential_1(x[0, 10], c[-0.1])" + // + // It defines en exponential exp(c * x). However, in HS3 the exponential + // is defined as exp(-c * x), to RooFit would export these dictionaries + // to the JSON: + // + // { + // "name": "exponential_1", // HS3 exponential_dist with transformed parameter + // "type": "exponential_dist", + // "x": "x", + // "c": "c_exponential_inverted" + // }, + // { + // "name": "c_exponential_inverted", // transformation function created on-the-fly on export + // "type": "generic_function", + // "expression": "-c" + // } + // + // On import, we can directly take the non-transformed parameter, which is + // we check for the suffix and optionally remove it from the requested + // name next: + + const std::string constParamName = p["c"].val(); + const bool isInverted = endsWith(constParamName, "_exponential_inverted"); + const std::string suffixToRemove = isInverted ? "_exponential_inverted" : ""; + RooAbsReal *c = tool->request(removeSuffix(constParamName, suffixToRemove), name); + + tool->wsEmplace(name, *x, *c, !isInverted); + + return true; + } +}; + +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 + { + std::string name(RooJSONFactoryWSTool::name(p)); + bool has_cov = p.has_child("covariances"); + bool has_corr = p.has_child("correlations") && p.has_child("standard_deviations"); + if (!has_cov && !has_corr) { + RooJSONFactoryWSTool::error("no covariances or correlations+standard_deviations given in '" + name + "'"); + } + + TMatrixDSym covmat; + + if (has_cov) { + int n = p["covariances"].num_children(); + int i = 0; + covmat.ResizeTo(n, n); + for (const auto &row : p["covariances"].children()) { + int j = 0; + for (const auto &val : row.children()) { + covmat(i, j) = val.val_double(); + ++j; + } + ++i; + } + } else { + std::vector variances; + for (const auto &v : p["standard_deviations"].children()) { + variances.push_back(v.val_double()); + } + covmat.ResizeTo(variances.size(), variances.size()); + int i = 0; + for (const auto &row : p["correlations"].children()) { + int j = 0; + for (const auto &val : row.children()) { + covmat(i, j) = val.val_double() * variances[i] * variances[j]; + ++j; + } + ++i; + } + } + tool->wsEmplace(name, tool->requestArgList(p, "x"), + tool->requestArgList(p, "mean"), covmat); + return true; + } +}; + +class ParamHistFuncFactory : public RooFit::JSONIO::Importer { +public: + bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override + { + std::string name(RooJSONFactoryWSTool::name(p)); + RooArgList varList = tool->requestArgList(p, "variables"); + if (!p.has_child("axes")) { + std::stringstream ss; + ss << "No axes given in '" << name << "'" + << ". Using default binning (uniform; nbins=100). If needed, export the Workspace to JSON with a newer " + << "Root version that supports custom ParamHistFunc binnings(>=6.38.00)." << std::endl; + RooJSONFactoryWSTool::warning(ss.str()); + tool->wsEmplace(name, varList, tool->requestArgList(p, "parameters")); + return true; + } + tool->wsEmplace(name, readBinning(p, varList), tool->requestArgList(p, "parameters")); + return true; + } + +private: + RooArgList readBinning(const JSONNode &topNode, const RooArgList &varList) const + { + // Temporary map from variable name → RooRealVar + std::map> varMap; + + // Build variables from JSON + for (const JSONNode &node : topNode["axes"].children()) { + const std::string name = node["name"].val(); + std::unique_ptr obs; + + if (node.has_child("edges")) { + std::vector edges; + for (const auto &bound : node["edges"].children()) { + edges.push_back(bound.val_double()); + } + obs = std::make_unique(name.c_str(), name.c_str(), edges.front(), edges.back()); + RooBinning bins(obs->getMin(), obs->getMax()); + for (auto b : edges) + bins.addBoundary(b); + obs->setBinning(bins); + } else { + obs = std::make_unique(name.c_str(), name.c_str(), node["min"].val_double(), + node["max"].val_double()); + obs->setBins(node["nbins"].val_int()); + } + + varMap[name] = std::move(obs); + } + + // 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)); + if (!refVar) + continue; + + auto it = varMap.find(refVar->GetName()); + if (it != varMap.end()) { + vars.addOwned(std::move(it->second)); // preserve ownership + } + } + return vars; + } +}; + +class RooSplineFactory : public RooFit::JSONIO::Importer { +public: + bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override + { + const std::string name(RooJSONFactoryWSTool::name(p)); + + // Mandatory fields + if (!p.has_child("x")) { + RooJSONFactoryWSTool::error("no x given in '" + name + "'"); + } + if (!p.has_child("x0") || !p.has_child("y0")) { + RooJSONFactoryWSTool::error("no x0/y0 given in '" + name + "'"); + } + + RooAbsReal *x = tool->requestArg(p, "x"); + + // Optional fields (defaults follow RooSpline ctor defaults) + std::string algo = p.has_child("interpolation") ? p["interpolation"].val() : "poly3"; + int order = 0; + if (algo == "poly3") + order = 3; + else if (algo == "poly5") + order = 5; + else { + RooJSONFactoryWSTool::error("unsupported algo '" + algo + "' for RooSpline in '" + name + + "': allowed are 'poly3' and 'poly5'"); + } + const bool logx = p.has_child("logx") ? p["logx"].val_bool() : false; + const bool logy = p.has_child("logy") ? p["logy"].val_bool() : false; + + // Read knots + std::vector x0; + std::vector y0; + x0.reserve(p["x0"].num_children()); + y0.reserve(p["y0"].num_children()); + + for (const auto &v : p["x0"].children()) + x0.push_back(v.val_double()); + for (const auto &v : p["y0"].children()) + y0.push_back(v.val_double()); + + if (x0.size() != y0.size()) { + RooJSONFactoryWSTool::error("x0/y0 size mismatch in '" + name + "': x0 has " + std::to_string(x0.size()) + + ", y0 has " + std::to_string(y0.size())); + } + if (x0.size() < 2) { + RooJSONFactoryWSTool::error("need at least 2 knots in '" + name + "'"); + } + + // Construct RooSpline(name,title, x, x0, y0, order, logx, logy) + tool->wsEmplace<::RooSpline>(name.c_str(), *x, std::span(x0.data(), x0.size()), + std::span(y0.data(), y0.size()), order, logx, logy); + + return true; + } +}; + +/////////////////////////////////////////////////////////////////////////////////////////////////////// +// specialized exporter implementations +/////////////////////////////////////////////////////////////////////////////////////////////////////// +template +class RooAddPdfStreamer : public RooFit::JSONIO::Exporter { +public: + std::string const &key() const override; + bool exportObject(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &elem) const override + { + const RooArg_t *pdf = static_cast(func); + elem["type"] << key(); + RooJSONFactoryWSTool::fillSeq(elem["summands"], pdf->pdfList()); + RooJSONFactoryWSTool::fillSeq(elem["coefficients"], pdf->coefList()); + elem["extended"] << (pdf->extendMode() != RooArg_t::CanNotBeExtended); + return true; + } +}; + +class RooRealSumPdfStreamer : public RooFit::JSONIO::Exporter { +public: + std::string const &key() const override; + bool exportObject(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &elem) const override + { + const RooRealSumPdf *pdf = static_cast(func); + elem["type"] << key(); + RooJSONFactoryWSTool::fillSeq(elem["samples"], pdf->funcList()); + RooJSONFactoryWSTool::fillSeq(elem["coefficients"], pdf->coefList()); + elem["extended"] << (pdf->extendMode() != RooAbsPdf::CanNotBeExtended); + return true; + } +}; + +class RooRealSumFuncStreamer : public RooFit::JSONIO::Exporter { +public: + std::string const &key() const override; + bool exportObject(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &elem) const override + { + const RooRealSumFunc *pdf = static_cast(func); + elem["type"] << key(); + RooJSONFactoryWSTool::fillSeq(elem["samples"], pdf->funcList()); + RooJSONFactoryWSTool::fillSeq(elem["coefficients"], pdf->coefList()); + return true; + } +}; + +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 { +public: + std::string const &key() const override; + bool exportObject(RooJSONFactoryWSTool *tool, const RooAbsArg *func, JSONNode &elem) const override + { + const RooHistPdf *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 RooHistPdfFactory : 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 RooBinSamplingPdfStreamer : public RooFit::JSONIO::Exporter { +public: + std::string const &key() const override; + bool exportObject(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &elem) const override + { + const RooBinSamplingPdf *pdf = static_cast(func); + elem["type"] << key(); + elem["pdf"] << pdf->pdf().GetName(); + elem["observable"] << pdf->observable().GetName(); + elem["epsilon"] << pdf->epsilon(); + return true; + } +}; + +class RooBinWidthFunctionStreamer : public RooFit::JSONIO::Exporter { +public: + std::string const &key() const override; + bool exportObject(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &elem) const override + { + const RooBinWidthFunction *pdf = static_cast(func); + elem["type"] << key(); + elem["histogram"] << pdf->histFunc().GetName(); + elem["divideByBinWidth"] << pdf->divideByBinWidth(); + return true; + } +}; + +template +class RooFormulaArgStreamer : public RooFit::JSONIO::Exporter { +public: + std::string const &key() const override; + bool exportObject(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &elem) const override + { + const RooArg_t *pdf = static_cast(func); + elem["type"] << key(); + TString expression(pdf->expression()); + cleanExpression(expression); + // If the tokens follow the "x[#]" convention, the square braces enclosing each number + // ensures that there is a unique mapping between the token and parameter name + // If the tokens follow the "@#" convention, the numbers are not enclosed by braces. + // So there may be tokens with numbers whose lower place value forms a subset string of ones with a higher place + // value, e.g. "@1" is a subset of "@10". So the names of these parameters must be applied descending from the + // highest place value in order to ensure each parameter name is uniquely applied to its token. + for (size_t idx = pdf->nParameters(); idx--;) { + const RooAbsArg *par = pdf->getParameter(idx); + expression.ReplaceAll(("x[" + std::to_string(idx) + "]").c_str(), par->GetName()); + expression.ReplaceAll(("@" + std::to_string(idx)).c_str(), par->GetName()); + } + elem["expression"] << expression.Data(); + return true; + } + +private: + void cleanExpression(TString &expr) const + { + expr.ReplaceAll("TMath::Exp", "exp"); + expr.ReplaceAll("TMath::Min", "min"); + expr.ReplaceAll("TMath::Max", "max"); + expr.ReplaceAll("TMath::Log", "log"); + expr.ReplaceAll("TMath::Cos", "cos"); + expr.ReplaceAll("TMath::Sin", "sin"); + expr.ReplaceAll("TMath::Sqrt", "sqrt"); + expr.ReplaceAll("TMath::Power", "pow"); + expr.ReplaceAll("TMath::Erf", "erf"); + } +}; +// Write the "x" reference and the coefficient list for polynomial-like +// pdfs/funcs, including the implicit defaults below "lowestOrder" so that the +// output is self-documenting. +template +void writePolynomialBody(const Pdf *pdf, JSONNode &elem) +{ + elem["x"] << pdf->x().GetName(); + auto &coefs = elem["coefficients"].set_seq(); + for (int i = 0; i < pdf->lowestOrder(); ++i) { + coefs.append_child() << (i == 0 ? "1.0" : "0.0"); + } + for (const auto &coef : pdf->coefList()) { + coefs.append_child() << coef->GetName(); + } +} + +template +class RooPolynomialStreamer : 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 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; + bool exportObject(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &elem) const override + { + auto *pdf = static_cast(func); + elem["type"] << key(); + elem["x"] << pdf->getX().GetName(); + elem["mean"] << pdf->getMean().GetName(); + elem["integer"] << !pdf->getNoRounding(); + return true; + } +}; + +class RooDecayStreamer : public RooFit::JSONIO::Exporter { +public: + std::string const &key() const override; + bool exportObject(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &elem) const override + { + auto *pdf = static_cast(func); + elem["type"] << key(); + elem["t"] << pdf->getT().GetName(); + elem["tau"] << pdf->getTau().GetName(); + elem["resolutionModel"] << pdf->getModel().GetName(); + elem["decayType"] << pdf->getDecayType(); + + return true; + } +}; + +class RooTruthModelStreamer : public RooFit::JSONIO::Exporter { +public: + std::string const &key() const override; + bool exportObject(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &elem) const override + { + auto *pdf = static_cast(func); + elem["type"] << key(); + elem["x"] << pdf->convVar().GetName(); + + return true; + } +}; + +class RooGaussModelStreamer : public RooFit::JSONIO::Exporter { +public: + std::string const &key() const override; + bool exportObject(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &elem) const override + { + auto *pdf = static_cast(func); + elem["type"] << key(); + elem["x"] << pdf->convVar().GetName(); + elem["mean"] << pdf->getMean().GetName(); + elem["sigma"] << pdf->getSigma().GetName(); + return true; + } +}; + +class RooLogNormalStreamer : public RooFit::JSONIO::Exporter { +public: + std::string const &key() const override; + bool exportObject(RooJSONFactoryWSTool *tool, const RooAbsArg *func, JSONNode &elem) const override + { + auto *pdf = static_cast(func); + + elem["type"] << key(); + elem["x"] << pdf->getX().GetName(); + + auto &m0 = pdf->getMedian(); + auto &k = pdf->getShapeK(); + + if (pdf->useStandardParametrization()) { + elem["mu"] << m0.GetName(); + elem["sigma"] << k.GetName(); + } else { + elem["mu"] << tool->exportTransformed(&m0, "_lognormal_log", "log(%s)"); + elem["sigma"] << tool->exportTransformed(&k, "_lognormal_log", "log(%s)"); + } + + return true; + } +}; + +class RooExponentialStreamer : public RooFit::JSONIO::Exporter { +public: + std::string const &key() const override; + bool exportObject(RooJSONFactoryWSTool *tool, const RooAbsArg *func, JSONNode &elem) const override + { + auto *pdf = static_cast(func); + elem["type"] << key(); + elem["x"] << pdf->variable().GetName(); + auto &c = pdf->coefficient(); + if (pdf->negateCoefficient()) { + elem["c"] << c.GetName(); + } else { + elem["c"] << tool->exportTransformed(&c, "_exponential_inverted", "-%s"); + } + + return true; + } +}; + +class RooMultiVarGaussianStreamer : public RooFit::JSONIO::Exporter { +public: + std::string const &key() const override; + bool exportObject(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &elem) const override + { + auto *pdf = static_cast(func); + elem["type"] << key(); + RooJSONFactoryWSTool::fillSeq(elem["x"], pdf->xVec()); + RooJSONFactoryWSTool::fillSeq(elem["mean"], pdf->muVec()); + elem["covariances"].fill_mat(pdf->covarianceMatrix()); + return true; + } +}; + +class RooTFnBindingStreamer : public RooFit::JSONIO::Exporter { +public: + std::string const &key() const override; + bool exportObject(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &elem) const override + { + auto *pdf = static_cast(func); + elem["type"] << key(); + + TString formula(pdf->function().GetExpFormula()); + formula.ReplaceAll("x", pdf->observables()[0].GetName()); + formula.ReplaceAll("y", pdf->observables()[1].GetName()); + formula.ReplaceAll("z", pdf->observables()[2].GetName()); + for (size_t i = 0; i < pdf->parameters().size(); ++i) { + TString pname(TString::Format("[%d]", (int)i)); + formula.ReplaceAll(pname, pdf->parameters()[i].GetName()); + } + elem["expression"] << formula.Data(); + return true; + } +}; + +class RooDerivativeStreamer : public RooFit::JSONIO::Exporter { +public: + std::string const &key() const override; + bool exportObject(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &elem) const override + { + auto *pdf = static_cast(func); + elem["type"] << key(); + elem["x"] << pdf->getX().GetName(); + elem["function"] << pdf->getFunc().GetName(); + if (!pdf->getNset().empty()) { + RooJSONFactoryWSTool::fillSeq(elem["normalization"], pdf->getNset()); + } + elem["order"] << pdf->order(); + elem["eps"] << pdf->eps(); + return true; + } +}; + +class RooRealIntegralStreamer : public RooFit::JSONIO::Exporter { +public: + std::string const &key() const override; + bool exportObject(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &elem) const override + { + auto *integral = static_cast(func); + elem["type"] << key(); + std::string integrand = integral->integrand().GetName(); + // elem["integrand"] << RooJSONFactoryWSTool::sanitizeName(integrand); + elem["integrand"] << integrand; + if (integral->intRange()) { + elem["domain"] << integral->intRange(); + } + RooJSONFactoryWSTool::fillSeq(elem["variables"], integral->intVars()); + if (RooArgSet const *funcNormSet = integral->funcNormSet()) { + RooJSONFactoryWSTool::fillSeq(elem["normalization"], *funcNormSet); + } + return true; + } +}; + +class RooFFTConvPdfStreamer : public RooFit::JSONIO::Exporter { +public: + std::string const &key() const override; + bool exportObject(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &elem) const override + { + auto *pdf = static_cast(func); + elem["type"] << key(); + if (auto convFunc = pdf->getPdfConvVar()) { + elem["conv_func"] << convFunc->GetName(); + } + elem["conv_var"] << pdf->getConvVar().GetName(); + elem["pdf1"] << pdf->getPdf1().GetName(); + elem["pdf2"] << pdf->getPdf2().GetName(); + elem["ipOrder"] << pdf->getInterpolationOrder(); + return true; + } +}; + +class RooExtendPdfStreamer : public RooFit::JSONIO::Exporter { +public: + std::string const &key() const override; + bool exportObject(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &elem) const override + { + auto *pdf = static_cast(func); + elem["type"] << key(); + if (auto rangeName = pdf->getRangeName()) { + elem["range"] << rangeName->GetName(); + } + elem["pdf"] << pdf->pdf().GetName(); + elem["norm"] << pdf->getN().GetName(); + return true; + } +}; + +class ParamHistFuncStreamer : public RooFit::JSONIO::Exporter { +public: + std::string const &key() const override; + bool exportObject(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &elem) const override + { + auto *pdf = static_cast(func); + elem["type"] << key(); + RooJSONFactoryWSTool::fillSeq(elem["variables"], pdf->dataVars()); + RooJSONFactoryWSTool::fillSeq(elem["parameters"], pdf->paramList()); + writeBinningInfo(pdf, elem); + return true; + } + +private: + void writeBinningInfo(const ParamHistFunc *pdf, JSONNode &elem) const + { + auto &observablesNode = elem["axes"].set_seq(); + // axes have to be ordered to get consistent bin indices + for (auto *var : static_range_cast(pdf->dataVars())) { + std::string name = var->GetName(); + RooJSONFactoryWSTool::testValidName(name, false); + JSONNode &obsNode = observablesNode.append_child().set_map(); + obsNode["name"] << name; + auto const &binning = var->getBinning(); + if (binning.isUniform()) { + obsNode["min"] << var->getMin(); + obsNode["max"] << var->getMax(); + obsNode["nbins"] << var->getBins(); + } else { + auto &edges = obsNode["edges"].set_seq(); + edges.append_child() << binning.binLow(0); + for (int i = 0; i < binning.numBins(); ++i) { + edges.append_child() << binning.binHigh(i); + } + } + } + } +}; + +class RooSplineStreamer : public RooFit::JSONIO::Exporter { +public: + std::string const &key() const override; + + bool exportObject(RooJSONFactoryWSTool *, const RooAbsArg *func, RooFit::Detail::JSONNode &elem) const override + { + auto const *rs = static_cast(func); + + elem["type"] << key(); + + // Independent variable + elem["x"] << rs->x().GetName(); + + // Spline configuration + // Canonical algo for RooSpline + elem["interpolation"] << (rs->order() == 5 ? "poly5" : "poly3"); + elem["logx"] << rs->logx(); + elem["logy"] << rs->logy(); + + // Serialize knots as primitive arrays + TSpline const &sp = rs->spline(); + auto &x0 = elem["x0"].set_seq(); + auto &y0 = elem["y0"].set_seq(); + + const int np = sp.GetNp(); + for (int i = 0; i < np; ++i) { + double xk = 0.0, yk = 0.0; + sp.GetKnot(i, xk, yk); + x0.append_child() << xk; + y0.append_child() << yk; + } + + return true; + } +}; + +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("misc", "ROOT_internal", "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["misc"]["ROOT_internal"]["selfNormalized"] << true; + + return true; + } +}; + +#define DEFINE_EXPORTER_KEY(class_name, name) \ + std::string const &class_name::key() const \ + { \ + const static std::string keystring = name; \ + return keystring; \ + } + +template <> +DEFINE_EXPORTER_KEY(RooAddPdfStreamer, "mixture_dist"); +template <> +DEFINE_EXPORTER_KEY(RooAddPdfStreamer, "mixture_model"); +DEFINE_EXPORTER_KEY(RooBinSamplingPdfStreamer, "binsampling"); +DEFINE_EXPORTER_KEY(RooWrapperPdfStreamer, "density_function_dist"); +DEFINE_EXPORTER_KEY(RooBinWidthFunctionStreamer, "binwidth"); +DEFINE_EXPORTER_KEY(RooLegacyExpPolyStreamer, "legacy_exp_poly_dist"); +DEFINE_EXPORTER_KEY(RooExponentialStreamer, "exponential_dist"); +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"); +DEFINE_EXPORTER_KEY(RooLogNormalStreamer, "lognormal_dist"); +DEFINE_EXPORTER_KEY(RooMultiVarGaussianStreamer, "multivariate_normal_dist"); +DEFINE_EXPORTER_KEY(RooPoissonStreamer, "poisson_dist"); +DEFINE_EXPORTER_KEY(RooDecayStreamer, "decay_dist"); +DEFINE_EXPORTER_KEY(RooTruthModelStreamer, "truth_model_function"); +DEFINE_EXPORTER_KEY(RooGaussModelStreamer, "gauss_model_function"); +template <> +DEFINE_EXPORTER_KEY(RooPolynomialStreamer, "polynomial_dist"); +template <> +DEFINE_EXPORTER_KEY(RooPolynomialStreamer, "polynomial"); +DEFINE_EXPORTER_KEY(RooRealSumFuncStreamer, "weighted_sum"); +DEFINE_EXPORTER_KEY(RooRealSumPdfStreamer, "weighted_sum_dist"); +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, "rate_extended_dist"); +DEFINE_EXPORTER_KEY(ParamHistFuncStreamer, "step"); +DEFINE_EXPORTER_KEY(RooSplineStreamer, "spline"); + +/////////////////////////////////////////////////////////////////////////////////////////////////////// +// instantiate all importers and exporters +/////////////////////////////////////////////////////////////////////////////////////////////////////// + +STATIC_EXECUTE([]() { + using namespace RooFit::JSONIO; + + registerImporter("density_function_dist"); + registerImporter("rate_extended_dist"); + registerImporter("product", false); + registerImporter("product_dist", false); + registerImporter("sum", false); + registerImporter("mixture_dist", false); + registerImporter("mixture_model", false); + registerImporter("binsampling_dist", false); + registerImporter("binwidth", 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("lognormal_dist", false); + registerImporter("multivariate_normal_dist", false); + registerImporter("poisson_dist", false); + registerImporter("decay_dist", false); + registerImporter("truth_model_function", false); + registerImporter("gauss_model_function", false); + registerImporter>("polynomial_dist", false); + registerImporter>("polynomial", false); + registerImporter("weighted_sum_dist", false); + registerImporter("weighted_sum", false); + registerImporter("integral", false); + registerImporter("derivative", false); + registerImporter("fft_conv_pdf", false); + registerImporter("extend_pdf", false); + 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(RooExponential::Class(), false); + registerExporter>(RooFormulaVar::Class(), false); + registerExporter>(RooGenericPdf::Class(), false); + registerExporter(RooHistFunc::Class(), false); + registerExporter(RooHistPdf::Class(), false); + registerExporter(RooLognormal::Class(), false); + registerExporter(RooMultiVarGaussian::Class(), false); + registerExporter(RooPoisson::Class(), false); + registerExporter(RooDecay::Class(), false); + registerExporter(RooTruthModel::Class(), false); + registerExporter(RooGaussModel::Class(), false); + registerExporter>(RooPolynomial::Class(), false); + registerExporter>(RooPolyVar::Class(), false); + registerExporter(RooRealSumFunc::Class(), false); + registerExporter(RooRealSumPdf::Class(), false); + registerExporter(RooTFnBinding::Class(), false); + registerExporter(RooRealIntegral::Class(), false); + registerExporter(RooDerivative::Class(), false); + registerExporter(RooFFTConvPdf::Class(), false); + registerExporter(RooExtendPdf::Class(), false); + registerExporter(ParamHistFunc::Class(), false); + registerExporter(RooSpline::Class(), false); +}); + +} // namespace diff --git a/roofit/hs3/src/JSONIOUtils.cxx b/roofit/hs3/src/JSONIOUtils.cxx index dd634f469444e..6e1687f5d2192 100644 --- a/roofit/hs3/src/JSONIOUtils.cxx +++ b/roofit/hs3/src/JSONIOUtils.cxx @@ -1,6 +1,9 @@ #include "JSONIOUtils.h" #include +#include +#include +#include #include #include @@ -167,3 +170,73 @@ std::vector readEdges(const JSONNode &node, std::string const &context) } return edges; } + +namespace { + +std::vector numBinsPerVariable(RooAbsCollection const &vars) +{ + std::vector nBins; + for (RooAbsArg const *var : vars) { + nBins.push_back(static_cast(dynamic_cast(*var).numBins())); + } + return nBins; +} + +// The position in TH1 order (first variable fastest) of the bin with position `k` in RooDataHist/HS3 order (last +// variable fastest). +std::size_t th1Position(std::size_t k, std::vector const &nBins) +{ + std::size_t th1 = 0; + std::size_t stride = 1; + std::vector indices(nBins.size()); + for (std::size_t j = nBins.size(); j-- > 0;) { + indices[j] = k % nBins[j]; + k /= nBins[j]; + } + for (std::size_t j = 0; j < nBins.size(); ++j) { + th1 += indices[j] * stride; + stride *= nBins[j]; + } + return th1; +} + +} // namespace + +std::vector paramHistFuncParamsInBinOrder(ParamHistFunc const &func) +{ + std::vector out; + for (int k = 0; k < func.numBins(); ++k) { + out.push_back(&func.getParameter(k)); + } + return out; +} + +RooArgList toParamHistFuncOrder(RooArgList const ¶msInBinOrder, RooAbsCollection const &vars) +{ + const std::vector nBins = numBinsPerVariable(vars); + std::vector th1Order(paramsInBinOrder.size(), nullptr); + for (std::size_t k = 0; k < paramsInBinOrder.size(); ++k) { + const std::size_t t = th1Position(k, nBins); + if (t >= th1Order.size()) { + RooJSONFactoryWSTool::error("the number of parameters does not match the number of bins"); + } + th1Order[t] = paramsInBinOrder.at(k); + } + RooArgList out; + for (RooAbsArg *arg : th1Order) { + out.add(*arg); + } + return out; +} + +void checkParamHistFuncOrderIsDefined(std::string const &hs3Version, RooAbsCollection const &vars, + std::string const &context) +{ + if (vars.size() >= 2 && isHS3VersionBefore(hs3Version, 0, 3)) { + RooJSONFactoryWSTool::error(context + + " has a multi-dimensional binning. Documents before HS3 0.3 don't define the order " + "of its per-bin parameters (RooFit wrote them in TH1 order, while the bin " + "contents use the order of HS3 0.3, last axis fastest). Please re-export the " + "workspace with a newer ROOT version."); + } +} diff --git a/roofit/hs3/src/JSONIOUtils.h b/roofit/hs3/src/JSONIOUtils.h index 2ac14662be7dc..8ea1e85213719 100644 --- a/roofit/hs3/src/JSONIOUtils.h +++ b/roofit/hs3/src/JSONIOUtils.h @@ -1,11 +1,45 @@ #ifndef JSONIOUtils_h #define JSONIOUtils_h +#include +#include +#include #include #include #include +class RooAbsArg; class RooAbsBinning; +class RooAbsCollection; +class RooAbsPdf; +class RooAbsReal; +class RooArgList; +class ParamHistFunc; + +// A RooProdPdf with this attribute represents an HS3 `array_dist`: a product of one-dimensional distributions of one +// type, each over its own observable. The factors are the components, imported as individual pdfs named after their +// observable. +inline constexpr const char *arrayDistTag = "hs3_array_dist"; + +// True if `arg` is a RooConstVar whose name is a plain number (i.e. a literal constant). +bool isLiteralConstVar(RooAbsArg const &arg); + +// The HS3 type and the observable of `pdf` as a component of an array distribution (RooPoisson, RooGaussian or +// RooLognormal in the standard parametrization, over a named observable), or std::nullopt if `pdf` cannot be a +// component. +// If `observable` names the mean of a RooGaussian, the Gaussian is written with x and mean swapped (see +// RooJSONFactoryWSTool::setConstraintObservable()). +struct ArrayComponentInfo { + std::string type; + std::string observable; +}; +std::optional arrayComponentInfo(RooAbsPdf const &pdf, std::string const *observable = nullptr); + +// Write the array distribution `name` with the given components into `elem`. `observableOf` gives the observable +// recorded for a component by name, if any. +void writeArrayDist(RooFit::Detail::JSONNode &elem, std::string const &name, + std::vector const &components, + std::map const &observableOf); bool startsWith(std::string_view str, std::string_view prefix); bool endsWith(std::string_view str, std::string_view suffix); @@ -17,4 +51,13 @@ bool parseHS3Version(std::string const &version, int &major, int &minor); bool isHS3VersionBefore(std::string const &version, int major, int minor); std::vector readEdges(const RooFit::Detail::JSONNode &node, std::string const &context); +// ParamHistFunc orders its parameters like the bins of a TH1 (first variable fastest), while HS3 lays out per-bin +// arrays with the last axis varying fastest, like RooDataHist. These helpers convert at the HS3 boundary. +std::vector paramHistFuncParamsInBinOrder(ParamHistFunc const &func); +RooArgList toParamHistFuncOrder(RooArgList const ¶msInBinOrder, RooAbsCollection const &vars); +// Documents before HS3 0.3 wrote the parameters of multi-dimensional ParamHistFuncs in TH1 order. No such documents +// are known; they are rejected instead of being guessed. +void checkParamHistFuncOrderIsDefined(std::string const &hs3Version, RooAbsCollection const &vars, + std::string const &context); + #endif diff --git a/roofit/hs3/src/RooFitHS3_wsexportkeys.cxx b/roofit/hs3/src/RooFitHS3_wsexportkeys.cxx index fec22e4350490..6e42794def03e 100644 --- a/roofit/hs3/src/RooFitHS3_wsexportkeys.cxx +++ b/roofit/hs3/src/RooFitHS3_wsexportkeys.cxx @@ -49,19 +49,6 @@ auto RooFitHS3_wsexportkeys = R"({ "coefList": "coefficients" } }, - "RooCrystalBall": { - "type": "crystalball_doublesided_dist", - "proxies": { - "alphaL": "alpha_L", - "alphaR": "alpha_R", - "nL": "n_L", - "nR": "n_R", - "x": "m", - "x0": "m0", - "sigmaL": "sigma_L", - "sigmaR": "sigma_R" - } - }, "RooEffProd": { "type": "efficiency_product_pdf_dist", "proxies": { diff --git a/roofit/hs3/src/RooFitHS3_wsfactoryexpressions.cxx b/roofit/hs3/src/RooFitHS3_wsfactoryexpressions.cxx index a938956e564c2..82fbf349da660 100644 --- a/roofit/hs3/src/RooFitHS3_wsfactoryexpressions.cxx +++ b/roofit/hs3/src/RooFitHS3_wsfactoryexpressions.cxx @@ -139,12 +139,6 @@ auto RooFitHS3_wsfactoryexpressions = R"({ "arguments": [ "x" ] - }, - "crystalball_doublesided_dist": { - "class": "RooCrystalBall", - "arguments": [ - "m", "m0", "sigma_L", "sigma_R", "alpha_L", "n_L", "alpha_R", "n_R" - ] } })"; diff --git a/roofit/hs3/src/RooJSONFactoryWSTool.cxx b/roofit/hs3/src/RooJSONFactoryWSTool.cxx index 9f372161b8e7c..826027bf7b099 100644 --- a/roofit/hs3/src/RooJSONFactoryWSTool.cxx +++ b/roofit/hs3/src/RooJSONFactoryWSTool.cxx @@ -20,9 +20,11 @@ #include #include #include +#include #include #include #include +#include #include #include #include @@ -34,6 +36,10 @@ #include #include +#include +#include +#include +#include #include #include #include @@ -106,7 +112,7 @@ tool.writedoc("hs3.tex") ~~~ */ -constexpr auto hs3VersionTag = "0.2"; +constexpr auto hs3VersionTag = "0.3"; using RooFit::Detail::JSONNode; using RooFit::Detail::JSONTree; @@ -367,8 +373,9 @@ std::vector> generateBinIndices(const RooArgSet &vars) std::vector> combinations; std::vector vars_numbins; vars_numbins.reserve(vars.size()); - for (const auto *absv : static_range_cast(vars)) { - vars_numbins.push_back(absv->getBins()); + for (const RooAbsArg *arg : vars) { + // Continuous variables and categories (one bin per state) + vars_numbins.push_back(dynamic_cast(*arg).numBins()); } std::vector curr_comb(vars.size()); ::genIndicesHelper(combinations, curr_comb, vars_numbins, 0); @@ -381,13 +388,6 @@ JSONNode const *findRooFitInternal(JSONNode const &node, Keys_t const &...keys) return node.find("misc", "ROOT_internal", keys...); } -// True if `arg` is a RooConstVar whose name is a plain number (i.e. a literal constant). -bool isLiteralConstVar(RooAbsArg const &arg) -{ - bool isRooConstVar = dynamic_cast(&arg); - return isRooConstVar && isNumber(arg.GetName()); -} - // Export the string attributes and tags of `arg` into the ROOT-internal attributes node. void exportAttributes(const RooAbsArg *arg, JSONNode &rootnode) { @@ -421,8 +421,8 @@ void exportAttributes(const RooAbsArg *arg, JSONNode &rootnode) } if (!arg->attributes().empty()) { for (auto const &attr : arg->attributes()) { - // Skip some RooFit internals - if (attr == "SnapShot_ExtRefClone" || attr == "RooRealConstant_Factory_Object") + // Skip some RooFit internals. The array distribution tag is restored from the type on import. + if (attr == "SnapShot_ExtRefClone" || attr == "RooRealConstant_Factory_Object" || attr == arrayDistTag) continue; initializeNode(); (*node)["tags"].set_seq().append_child() << attr; @@ -435,7 +435,27 @@ void getObservables(RooWorkspace const &ws, const JSONNode &node, RooAbsCollecti { for (const auto &p : node["axes"].children()) { std::string name(RooJSONFactoryWSTool::name(p)); - if (ws.var(name)) { + if (p.has_child("categories")) { + RooCategory *cat = ws.cat(name); + if (!cat) { + RooJSONFactoryWSTool::error("The categorical observable \"" + name + + "\" could not be found in the workspace! Categorical observables are " + "declared by a categorical axis in a domain."); + } + std::set labels; + for (auto const &label : p["categories"].children()) { + labels.insert(label.val()); + } + std::set states; + for (auto const &state : *cat) { + states.insert(state.first); + } + if (labels != states) { + RooJSONFactoryWSTool::error("The categories of the axis \"" + name + + "\" differ from the states of the category \"" + name + "\""); + } + out.add(*cat); + } else if (ws.var(name)) { out.add(*ws.var(name)); } else { RooJSONFactoryWSTool::error("The observable \"" + name + "\" could not be found in the workspace!"); @@ -470,6 +490,40 @@ void applyUnbinnedAxisRanges(const JSONNode &p, RooWorkspace &workspace, RooFit: } } +// Give the categories among `vars` the states of the categories with the same name in the workspace: the same +// label-to-index mapping and the same order. RooDataHist looks up the bin of a category by the order of the states of +// the category it is given, so data and workspace category have to agree on it. The workspace categories come from +// the domains; the labels must be the same, but the data may list them in a different order (HS3 identifies +// categories by name). readBinnedData() then permutes the bin contents accordingly. +void alignCategoriesWithWorkspace(RooArgSet &vars, RooWorkspace const &workspace, std::string const &dataName) +{ + for (RooAbsArg *arg : vars) { + auto *cat = dynamic_cast(arg); + if (!cat) { + continue; + } + RooCategory const *wsCat = workspace.cat(cat->GetName()); + if (!wsCat) { + continue; + } + std::vector labels; + for (int i = 0; i < cat->numTypes(); ++i) { + labels.push_back(cat->lookupName(i)); + } + if (static_cast(labels.size()) != wsCat->numTypes() || + std::any_of(labels.begin(), labels.end(), [&](auto const &l) { return !wsCat->hasLabel(l); })) { + RooJSONFactoryWSTool::error("the categories of the axis '" + std::string(cat->GetName()) + + "' of the dataset '" + dataName + "' differ from those of the category '" + + cat->GetName() + "'"); + } + cat->clear(); + for (unsigned int i = 0; i < wsCat->size(); ++i) { + auto const &state = wsCat->getOrdinal(i); + cat->defineType(state.first, state.second); + } + } +} + // Create a RooAbsData (binned or unbinned) from the JSON node `p`. std::unique_ptr loadData(const JSONNode &p, RooWorkspace &workspace, RooFit::JSONIO::Detail::Domains const *domains) @@ -481,7 +535,9 @@ loadData(const JSONNode &p, RooWorkspace &workspace, RooFit::JSONIO::Detail::Dom std::string const &type = p["type"].val(); if (type == "binned") { // binned - return RooJSONFactoryWSTool::readBinnedData(p, name, RooJSONFactoryWSTool::readAxes(p)); + RooArgSet axes = RooJSONFactoryWSTool::readAxes(p); + alignCategoriesWithWorkspace(axes, workspace, name); + return RooJSONFactoryWSTool::readBinnedData(p, name, axes); } else if (type == "unbinned") { // unbinned RooArgList varlist; @@ -513,8 +569,17 @@ loadData(const JSONNode &p, RooWorkspace &workspace, RooFit::JSONIO::Detail::Dom } std::size_t j = 0; for (auto const &pointj : point.children()) { - auto *v = static_cast(varlist.at(j)); - v->setVal(pointj.val_double()); + if (auto *cat = dynamic_cast(varlist.at(j))) { + // The entries of a categorical axis are the labels of the categories. + if (pointj.is_number() || !cat->hasLabel(pointj.val())) { + RooJSONFactoryWSTool::error("entry " + std::to_string(i) + " of the dataset '" + name + + "' has the value '" + pointj.val() + "' on the categorical axis '" + + cat->GetName() + "', which is not one of its categories"); + } + cat->setLabel(pointj.val()); + } else { + static_cast(varlist.at(j))->setVal(pointj.val_double()); + } ++j; } if (weightVals.size() > 0) { @@ -533,7 +598,7 @@ loadData(const JSONNode &p, RooWorkspace &workspace, RooFit::JSONIO::Detail::Dom // Import an analysis (likelihood + domains) as one or more ModelConfig objects into the workspace. void importAnalysis(const JSONNode &rootnode, const JSONNode &analysisNode, const JSONNode &likelihoodsNode, const JSONNode &domainsNode, RooWorkspace &workspace, - const std::vector> &datasets) + std::vector> &datasets) { // if this is a toplevel pdf, also create a modelConfig for it std::string const &analysisName = RooJSONFactoryWSTool::name(analysisNode); @@ -558,24 +623,99 @@ void importAnalysis(const JSONNode &rootnode, const JSONNode &analysisNode, cons if (!nllNode->has_child("data")) { throw std::runtime_error("likelihood node has no data attached!"); } - std::vector nllDistNames = valsToStringVec((*nllNode)["distributions"]); + std::vector const listedDistNames = valsToStringVec((*nllNode)["distributions"]); + std::vector const listedDataNames = valsToStringVec((*nllNode)["data"]); + + // The constraint terms of the HistFactory distributions of the likelihood enter through those distributions. + // Listed in "aux_distributions", they are not added again; paired with a data set in "distributions", they are + // not channels, and the data set gives the observed values of their observables. + std::vector> trees; + for (std::string const &name : listedDistNames) { + RooAbsPdf *listed = workspace.pdf(name); + trees.push_back(listed ? std::unique_ptr{listed->getComponents()} : nullptr); + } + auto isConstraintTermOfOthers = [&](RooAbsArg const &term, std::size_t self) { + auto isPartOf = [&term](RooArgSet const &tree) { + if (tree.find(term.GetName())) { + return true; + } + if (term.getAttribute(arrayDistTag)) { + if (auto const *array = dynamic_cast(&term)) { + for (RooAbsArg const *component : array->pdfList()) { + if (tree.find(component->GetName())) { + return true; + } + } + } + } + return false; + }; + for (std::size_t i = 0; i < trees.size(); ++i) { + if (i != self && trees[i] && isPartOf(*trees[i])) { + return true; + } + } + return false; + }; + RooArgSet extConstraints; for (auto &nameNode : (*nllNode)["aux_distributions"].children()) { if (RooAbsArg *extConstraint = workspace.arg(nameNode.val())) { - extConstraints.add(*extConstraint); + if (!isConstraintTermOfOthers(*extConstraint, listedDistNames.size())) { + extConstraints.add(*extConstraint); + } + } + } + + std::vector nllDistNames; + std::vector nllDataNames; + RooArgSet constraintObservables; + for (std::size_t i = 0; i < listedDistNames.size(); ++i) { + std::string const &distName = listedDistNames[i]; + std::string const dataName = i < listedDataNames.size() ? listedDataNames[i] : "0"; + RooAbsPdf *listed = workspace.pdf(distName); + if (!listed || !isConstraintTermOfOthers(*listed, i)) { + nllDistNames.push_back(distName); + nllDataNames.push_back(dataName); + continue; + } + for (auto &d : datasets) { + if (!d || d->GetName() != dataName) { + continue; + } + if (d->numEntries() != 1) { + throw std::runtime_error("the data set '" + dataName + "' paired with the constraint term '" + distName + + "' must have exactly one entry"); + } + for (RooAbsArg const *obs : *d->get(0)) { + auto *var = workspace.var(obs->GetName()); + auto const *value = dynamic_cast(obs); + if (!var || !value) { + RooJSONFactoryWSTool::warning("RooFitHS3: cannot set the observed value of '" + + std::string(obs->GetName()) + "' from the data set '" + dataName + + "': it is not a variable of the workspace"); + continue; + } + var->setVal(value->getVal()); + var->setConstant(true); + constraintObservables.add(*var); + } + // The values live in the variables now; the data set is not imported. + d.reset(); } } + RooArgSet observables; - for (auto &nameNode : (*nllNode)["data"].children()) { + for (std::size_t i = 0; i < nllDataNames.size(); ++i) { bool found = false; for (const auto &d : datasets) { - if (d->GetName() == nameNode.val()) { + if (d && d->GetName() == nllDataNames[i]) { found = true; observables.add(*d->get(), true); } } - if (nameNode.val() != "0" && !found) - throw std::runtime_error("dataset '" + nameNode.val() + "' cannot be found!"); + if (nllDataNames[i] != "0" && !found) + throw std::runtime_error("dataset '" + nllDataNames[i] + "' cannot be found!"); } JSONNode const *pdfNameNode = mcAuxNode ? mcAuxNode->find("pdfName") : nullptr; @@ -609,15 +749,35 @@ void importAnalysis(const JSONNode &rootnode, const JSONNode &analysisNode, cons if (!extConstraints.empty()) mc->SetExternalConstraints(extConstraints); - auto readArgSet = [&](std::string const &name) { - RooArgSet out; - for (auto const &child : analysisNode[name].children()) { - out.add(*workspace.arg(child.val())); + // HS3 v0.2.9 used the singular keys "parameter_of_interest" (a list) and "domain" (a single name). They are still + // accepted, but the plural keys take precedence. + auto findAnalysisKey = [&](const char *key, const char *legacyKey) -> JSONNode const * { + if (auto node = analysisNode.find(key)) { + return node; + } + if (auto node = analysisNode.find(legacyKey)) { + RooJSONFactoryWSTool::warning(std::string("RooFitHS3: the key \"") + legacyKey + "\" of analysis '" + + analysisName + "' is deprecated, please use \"" + key + "\" instead."); + return node; } - return out; + return nullptr; }; + JSONNode const *poiNode = findAnalysisKey("parameters_of_interest", "parameter_of_interest"); + JSONNode const *domainsNodeOfAnalysis = findAnalysisKey("domains", "domain"); + std::vector analysisDomains; + if (domainsNodeOfAnalysis && domainsNodeOfAnalysis->is_seq()) { + analysisDomains = valsToStringVec(*domainsNodeOfAnalysis); + } else if (domainsNodeOfAnalysis) { + analysisDomains.push_back(domainsNodeOfAnalysis->val()); + } - mc->SetParametersOfInterest(readArgSet("parameters_of_interest")); + RooArgSet pois; + if (poiNode) { + for (auto const &child : poiNode->children()) { + pois.add(*workspace.arg(child.val())); + } + } + mc->SetParametersOfInterest(pois); mc->SetObservables(observables); RooArgSet pars; pdf->getParameters(&observables, pars); @@ -629,8 +789,8 @@ void importAnalysis(const JSONNode &rootnode, const JSONNode &analysisNode, cons pdf->getAllConstraints(observables, mainPars, /*stripDisconnected*/ true); RooArgSet domainPars; - for (auto &domain : analysisNode["domains"].children()) { - const auto &thisDomain = RooJSONFactoryWSTool::findNamedChild(domainsNode, domain.val()); + for (std::string const &domain : analysisDomains) { + const auto &thisDomain = RooJSONFactoryWSTool::findNamedChild(domainsNode, domain); if (!thisDomain || !thisDomain->has_child("axes")) continue; for (auto &var : (*thisDomain)["axes"].children()) { @@ -652,6 +812,7 @@ void importAnalysis(const JSONNode &rootnode, const JSONNode &analysisNode, cons } } + globs.add(constraintObservables, /*silent=*/true); mc->SetGlobalObservables(globs); mc->SetNuisanceParameters(nps); @@ -686,6 +847,39 @@ void importAnalysis(const JSONNode &rootnode, const JSONNode &analysisNode, cons } } +// HistFactory modifiers pick the components that constrain their parameters out of an array distribution. By the +// standard the whole array enters the likelihood, so the components of a partially used array that nothing uses are +// reported. +void warnAboutUnusedArrayComponents(RooWorkspace const &ws) +{ + for (RooAbsArg *arg : ws.allPdfs()) { + auto const *array = dynamic_cast(arg); + if (!array || !array->getAttribute(arrayDistTag) || array->hasClients()) { + continue; + } + bool used = false; + std::string unused; + for (RooAbsArg const *component : array->pdfList()) { + bool usedElsewhere = false; + for (RooAbsArg const *client : component->clients()) { + if (client != array) { + usedElsewhere = true; + } + } + if (usedElsewhere) { + used = true; + } else { + unused += (unused.empty() ? "'" : ", '") + std::string(component->GetName()) + "'"; + } + } + if (used && !unused.empty()) { + RooJSONFactoryWSTool::warning("RooFitHS3: the components " + unused + " of the array distribution '" + + array->GetName() + + "' are not used by any distribution and are not part of any likelihood"); + } + } +} + void combinePdfs(const JSONNode &rootnode, RooWorkspace &ws) { auto *combinedPdfInfoNode = findRooFitInternal(rootnode, "combined_distributions"); @@ -825,8 +1019,10 @@ void checkHS3Version(std::string const &version) } } -/// HS3 requires the names of all top-level objects to be unique within a document. Throw if a name is used twice, or -/// only warn if not `strict`: RooFit wrote such documents before HS3 v0.3, e.g. for a range and a pdf of one name. +/// HS3 requires the names of all top-level objects to be unique within a document. If `strict`, throw if a name is +/// used twice within one section or by a distribution and a function. Other clashes, and all clashes if not +/// `strict`, only warn: RooFit's own export still produces them, e.g. for a range and a pdf of one name, and it +/// must be able to read its own output. void checkUniqueNames(JSONNode const &rootnode, bool strict) { std::map sectionOfName; @@ -850,7 +1046,11 @@ void checkUniqueNames(JSONNode const &rootnode, bool strict) // compatibility: analyses are not referenced by name, and likelihoods only from analyses, so there is no // ambiguity. auto isLegacyNamed = [](std::string const &s) { return s == "likelihoods" || s == "analyses"; }; - if (!strict || isLegacyNamed(firstSection) || isLegacyNamed(section)) { + // Within one section, and between distributions and functions (both can be referenced by name), a clash is + // ambiguous. Other clashes are produced by RooFit itself (a range, snapshot or data set named like a pdf). + auto isReferable = [](std::string const &s) { return s == "distributions" || s == "functions"; }; + const bool ambiguous = firstSection == section || (isReferable(firstSection) && isReferable(section)); + if (!strict || !ambiguous || isLegacyNamed(firstSection) || isLegacyNamed(section)) { RooJSONFactoryWSTool::warning("RooFitHS3: the name '" + name + "' is used in \"" + firstSection + "\" and in \"" + section + "\", which is not allowed by the HS3 standard."); continue; @@ -863,6 +1063,12 @@ void checkUniqueNames(JSONNode const &rootnode, bool strict) } // namespace +bool isLiteralConstVar(RooAbsArg const &arg) +{ + bool isRooConstVar = dynamic_cast(&arg); + return isRooConstVar && isNumber(arg.GetName()); +} + RooJSONFactoryWSTool::RooJSONFactoryWSTool(RooWorkspace &ws) : _workspace{ws} {} RooJSONFactoryWSTool::~RooJSONFactoryWSTool() {} @@ -1078,21 +1284,31 @@ std::string RooJSONFactoryWSTool::exportTransformed(const RooAbsReal *original, } // Export a single object `func` (pdf, function, variable or category) to the output JSON, recording its name in -// `exportedObjectNames` to avoid exporting it twice. -void RooJSONFactoryWSTool::exportObject(RooAbsArg const &func, std::set &exportedObjectNames) +// `_exportedObjectNames` to avoid exporting it twice. +void RooJSONFactoryWSTool::exportObject(RooAbsArg const &func) { std::string name = func.GetName(); // if this element was already exported, skip - if (exportedObjectNames.find(name) != exportedObjectNames.end()) + if (_exportedObjectNames.find(name) != _exportedObjectNames.end()) return; - exportedObjectNames.insert(name); + _exportedObjectNames.insert(name); + + // Exporters that group objects per top-level pdf (the HistFactory constraints) ask for the root of the recursion. + if (_exportDepth == 0) { + _exportRootName = name; + } + struct DepthGuard { + int &depth; + DepthGuard(int &d) : depth{d} { ++depth; } + ~DepthGuard() { --depth; } + } depthGuard{_exportDepth}; if (auto simPdf = dynamic_cast(&func)) { // RooSimultaneous is not used in the HS3 standard, we only export the // dependents and some ROOT internal information. - exportObjects(func.servers(), exportedObjectNames); + exportObjects(func.servers()); std::vector channelNames; for (auto const &item : simPdf->indexCat()) { @@ -1115,6 +1331,22 @@ void RooJSONFactoryWSTool::exportObject(RooAbsArg const &func, std::set(&func) || dynamic_cast(&func)) { exportVariable(&func, *_varsNode, true, false); return; + } else if (config().exportArrayDist && func.getAttribute(arrayDistTag)) { + // An array distribution (HS3 `array_dist`): its components are registered here, and the array is written at + // the end of the export, see exportAsArrayComponent(). + auto const *prod = dynamic_cast(&func); + if (!prod) { + error("'" + name + "' has the attribute '" + arrayDistTag + + "' of an array distribution, but it is not a RooProdPdf"); + } + for (RooAbsPdf *pdf : static_range_cast(prod->pdfList())) { + const std::string writtenIn = exportAsArrayComponent(*pdf, name); + if (writtenIn != name) { + error("cannot write '" + name + "' as an array distribution: its component '" + pdf->GetName() + + "' is already written as '" + writtenIn + "'"); + } + } + return; } auto &collectionNode = (*_rootnodeOutput)[dynamic_cast(&func) ? "distributions" : "functions"]; @@ -1142,9 +1374,9 @@ void RooJSONFactoryWSTool::exportObject(RooAbsArg const &func, std::setautoExportDependants()) { - exportObjects(func.servers(), exportedObjectNames); + exportObjects(func.servers()); } else { - exportObjects(_serversToExport, exportedObjectNames); + exportObjects(_serversToExport); } for (auto &s : _serversToDelete) { delete s; @@ -1216,10 +1448,137 @@ void RooJSONFactoryWSTool::exportObject(RooAbsArg const &func, std::setexportObject(*s, exportedObjectNames); + this->exportObject(*s); + } +} + +std::string RooJSONFactoryWSTool::exportAsArrayComponent(RooAbsPdf const &pdf, std::string const &arrayName) +{ + const std::string pdfName = pdf.GetName(); + if (auto it = _arrayOfComponent.find(pdfName); it != _arrayOfComponent.end()) { + return it->second; + } + // Already written on its own: it is referenced by its name. + if (_exportedObjectNames.find(pdfName) != _exportedObjectNames.end()) { + return pdfName; + } + if (_arrayComponents.find(arrayName) == _arrayComponents.end()) { + RooAbsArg const *existing = _workspace.arg(arrayName); + if (existing && !existing->getAttribute(arrayDistTag)) { + error("cannot write the array distribution '" + arrayName + + "': the workspace has another object of that name"); + } + _exportedObjectNames.insert(arrayName); + } + _arrayComponents[arrayName].push_back(&pdf); + _arrayOfComponent[pdfName] = arrayName; + _exportedObjectNames.insert(pdfName); + return arrayName; +} + +// Write the array distributions registered with exportAsArrayComponent(), in alphabetical order, each followed by the +// servers of its components. +void RooJSONFactoryWSTool::exportArrayDists() +{ + std::set written; + for (;;) { + auto next = std::find_if(_arrayComponents.begin(), _arrayComponents.end(), + [&](auto const &item) { return written.find(item.first) == written.end(); }); + if (next == _arrayComponents.end()) { + break; + } + written.insert(next->first); + // Exporting the servers may add to the registry, so the components are copied. + const std::vector components = next->second; + writeArrayDist(appendNamedChild((*_rootnodeOutput)["distributions"], next->first), next->first, components, + _constraintObservables); + for (RooAbsPdf const *pdf : components) { + exportObjects(pdf->servers()); + } + } +} + +void RooJSONFactoryWSTool::setConstraintObservable(RooAbsPdf const &pdf, std::string const &observable) +{ + _constraintObservables[pdf.GetName()] = observable; +} + +std::string const *RooJSONFactoryWSTool::constraintObservable(std::string const &pdfName) const +{ + auto found = _constraintObservables.find(pdfName); + return found == _constraintObservables.end() ? nullptr : &found->second; +} + +void RooJSONFactoryWSTool::recordConstraintTerms(std::string const &pdfName, std::vector const &terms) +{ + auto &out = _constraintTermsOfPdf[pdfName]; + for (auto const &term : terms) { + if (std::find(out.begin(), out.end(), term) == out.end()) { + out.push_back(term); + } } } +// Write the data set "_data" with the observed values of the observables of the constraint term `term` (an +// array distribution written in this export, or a single constraint pdf) and return true, or return false if the +// observables of the term are not all global observables of the model. +bool RooJSONFactoryWSTool::exportConstraintData(std::string const &term, RooArgSet const *globalObservables, + JSONNode &rootnode) +{ + std::vector components; + if (auto found = _arrayComponents.find(term); found != _arrayComponents.end()) { + components = found->second; + } else if (RooAbsPdf *pdf = _workspace.pdf(term)) { + components.push_back(pdf); + } + std::vector observables; + for (RooAbsPdf const *pdf : components) { + auto info = arrayComponentInfo(*pdf, constraintObservable(pdf->GetName())); + if (!info || !globalObservables || !globalObservables->find(info->observable.c_str())) { + return false; + } + observables.push_back(info->observable); + } + if (observables.empty()) { + return false; + } + const std::string dataName = term + "_data"; + if (!_writtenConstraintData.insert(dataName).second) { + return true; + } + if (_workspace.data(dataName)) { + error("cannot write the observed values of the constraint term '" + term + "' as the data set '" + dataName + + "': the workspace has a data set of that name"); + } + JSONNode &node = appendNamedChild(rootnode["data"], dataName); + node["type"] << "unbinned"; + auto &axes = node["axes"].set_seq(); + auto &entry = node["entries"].set_seq().append_child().set_seq(); + for (std::string const &name : observables) { + auto const *var = dynamic_cast(_workspace.arg(name)); + if (!var) { + error("the observable '" + name + "' of the constraint term '" + term + "' is not in the workspace"); + } + auto &axis = axes.append_child().set_map(); + axis["name"] << name; + // A range without a bound gets the largest representable value, so that other engines see a finite axis. + double min = -std::numeric_limits::max(); + double max = std::numeric_limits::max(); + if (auto const *rrv = dynamic_cast(var)) { + if (std::isfinite(rrv->getMin())) { + min = rrv->getMin(); + } + if (std::isfinite(rrv->getMax())) { + max = rrv->getMax(); + } + } + axis["min"] << min; + axis["max"] << max; + entry.append_child() << var->getVal(); + } + return true; +} + /** * @brief Import a function from the JSONNode into the workspace. * @@ -1354,6 +1713,11 @@ void RooJSONFactoryWSTool::exportAxis(JSONNode &obsNode, RooRealVar const &var) obsNode["name"] << name; writeAxisBinning(obsNode, var.getBinning()); + // Binned axes always carry their range, also if it is redundant with the "edges" of an irregular binning. + if (!obsNode.has_child("min")) { + obsNode["min"] << var.getBinning().lowBound(); + obsNode["max"] << var.getBinning().highBound(); + } } /** @@ -1591,7 +1955,28 @@ RooArgSet RooJSONFactoryWSTool::readAxes(const JSONNode &topNode) RooArgSet vars; for (JSONNode const &node : topNode["axes"].children()) { - if (node.has_child("edges")) { + if (node.has_child("categories")) { + // One state per category, in the order of the document, with the position as index. The bins of a + // RooDataHist then follow the order of "categories". + const std::string axisName = node["name"].val(); + for (const char *key : {"min", "max", "nbins", "edges"}) { + if (node.has_child(key)) { + error("the axis '" + axisName + "' is categorical and must not have \"" + key + "\""); + } + } + auto cat = std::make_unique(axisName.c_str(), axisName.c_str()); + int index = 0; + for (JSONNode const &label : node["categories"].children()) { + if (label.val().empty() || cat->hasLabel(label.val())) { + error("the categories of the axis '" + axisName + "' must be unique non-empty strings"); + } + cat->defineType(label.val(), index++); + } + if (index == 0) { + error("the categorical axis '" + axisName + "' has no categories"); + } + vars.addOwned(std::move(cat)); + } else if (node.has_child("edges")) { std::vector edges = readEdges(node["edges"], "axis '" + node["name"].val() + "'"); auto obs = std::make_unique(node["name"].val().c_str(), node["name"].val().c_str(), edges[0], edges[edges.size() - 1]); @@ -1612,6 +1997,82 @@ RooArgSet RooJSONFactoryWSTool::readAxes(const JSONNode &topNode) return vars; } +namespace { + +// For binned data in `n` with categorical axes, the document lays out the bins in the order of its "categories" (last +// axis fastest), while RooFit lays them out in the order of the states of the categories in `vars`. Returns, for each +// bin in RooFit's order, the position of that bin in the document. Empty if the orders agree. +std::vector binOrderOfDocument(const JSONNode &n, RooArgSet const &vars) +{ + JSONNode const *axes = n.find("axes"); + if (!axes || static_cast(axes->num_children()) != vars.size()) { + return {}; + } + std::vector nBins; + // For each axis, the position in the document of the bin with a given position in RooFit's order. + std::vector> axisPositions; + bool reordered = false; + std::size_t iAxis = 0; + for (JSONNode const &axis : axes->children()) { + RooAbsArg const *var = vars[iAxis++]; + auto const *cat = dynamic_cast(var); + const auto numBins = static_cast(dynamic_cast(*var).numBins()); + nBins.push_back(numBins); + std::vector positions(numBins); + for (std::size_t i = 0; i < numBins; ++i) { + positions[i] = i; + } + if (cat && axis.has_child("categories")) { + std::vector labels; + for (JSONNode const &label : axis["categories"].children()) { + labels.push_back(label.val()); + } + for (std::size_t i = 0; i < numBins; ++i) { + auto found = std::find(labels.begin(), labels.end(), cat->getOrdinal(i).first); + positions[i] = std::distance(labels.begin(), found); + reordered |= positions[i] != i; + } + } + axisPositions.push_back(std::move(positions)); + } + if (!reordered) { + return {}; + } + std::size_t total = 1; + for (std::size_t nb : nBins) { + total *= nb; + } + std::vector order(total); + for (std::size_t k = 0; k < total; ++k) { + // Decompose RooFit's bin index (last axis fastest) and recompose it with the document's positions. + std::size_t rest = k; + std::size_t docIndex = 0; + std::size_t stride = 1; + for (std::size_t j = nBins.size(); j-- > 0;) { + const std::size_t i = rest % nBins[j]; + rest /= nBins[j]; + docIndex += axisPositions[j][i] * stride; + stride *= nBins[j]; + } + order[k] = docIndex; + } + return order; +} + +std::vector permute(std::vector const &values, std::vector const &order) +{ + if (order.empty() || values.size() != order.size()) { + return values; + } + std::vector out(values.size()); + for (std::size_t k = 0; k < order.size(); ++k) { + out[k] = values[order[k]]; + } + return out; +} + +} // namespace + /** * @brief Read binned data from the JSONNode and create a RooDataHist object. * @@ -1652,17 +2113,20 @@ RooJSONFactoryWSTool::readBinnedData(const JSONNode &n, const std::string &name, ", bins=" + std::to_string(bins.size())); } auto dh = std::make_unique(name, name, vars); + const std::vector binOrder = binOrderOfDocument(n, vars); std::vector contentVals; contentVals.reserve(contents.num_children()); for (auto const &cont : contents.children()) { contentVals.push_back(cont.val_double()); } + contentVals = permute(contentVals, binOrder); std::vector errorVals; if (errors) { errorVals.reserve(errors->num_children()); for (auto const &err : errors->children()) { errorVals.push_back(err.val_double()); } + errorVals = permute(errorVals, binOrder); } // RooDataHist only starts to store squared weights when it sees a positive error, and then initializes them from // the weights so far. Trigger that before setting any weight, so that errors of zero stay zero. @@ -1793,13 +2257,52 @@ void RooJSONFactoryWSTool::exportSingleModelConfig(JSONNode &rootnode, RooFit::M nllNode["data"].append_child() << 0; } + // The additional terms of the likelihood: the external constraints of the ModelConfig, and the constraint terms + // of the HistFactory channels of the model. The latter are listed in "aux_distributions" with the observed values + // of their observables as constant parameters, or (config().exportConstraintData) in "distributions", each paired + // with a data set that holds the observed values. + std::vector auxTerms; + auto addAuxTerm = [&](std::string const &term) { + if (std::find(auxTerms.begin(), auxTerms.end(), term) == auxTerms.end()) { + auxTerms.push_back(term); + } + }; if (mc.GetExternalConstraints()) { - auto &extConstrNode = nllNode["aux_distributions"]; - extConstrNode.set_seq(); for (const auto &constr : *mc.GetExternalConstraints()) { - extConstrNode.append_child() << constr->GetName(); + addAuxTerm(constr->GetName()); + } + } + std::vector channels; + if (auto simPdf = dynamic_cast(pdf)) { + for (auto const &item : simPdf->indexCat()) { + channels.push_back(simPdf->getPdf(item.first)->GetName()); + } + } else { + channels.push_back(pdf->GetName()); + } + std::vector constraintTerms; + for (std::string const &channel : channels) { + auto found = _constraintTermsOfPdf.find(channel); + if (found == _constraintTermsOfPdf.end()) { + continue; + } + for (std::string const &term : found->second) { + if (std::find(constraintTerms.begin(), constraintTerms.end(), term) == constraintTerms.end()) { + constraintTerms.push_back(term); + } + } + } + for (std::string const &term : constraintTerms) { + if (config().exportConstraintData && exportConstraintData(term, mc.GetGlobalObservables(), rootnode)) { + nllNode["distributions"].append_child() << term; + nllNode["data"].append_child() << term + "_data"; + } else { + addAuxTerm(term); } } + if (!auxTerms.empty()) { + nllNode["aux_distributions"].fill_seq(auxTerms); + } auto writeList = [&](const char *name, RooArgSet const *args) { if (!args || !args->size()) @@ -1861,8 +2364,13 @@ void RooJSONFactoryWSTool::exportAllObjects(JSONNode &n) } } sortByName(allpdfs); - std::set exportedObjectNames; - exportObjects(allpdfs, exportedObjectNames); + _exportedObjectNames.clear(); + _arrayComponents.clear(); + _arrayOfComponent.clear(); + _constraintTermsOfPdf.clear(); + _writtenConstraintData.clear(); + _constraintObservables.clear(); + exportObjects(allpdfs); // export all toplevel functions std::vector allfuncs; @@ -1874,7 +2382,10 @@ void RooJSONFactoryWSTool::exportAllObjects(JSONNode &n) } } sortByName(allfuncs); - exportObjects(allfuncs, exportedObjectNames); + exportObjects(allfuncs); + + // The array distributions are written after the objects that use them, see exportAsArrayComponent(). + exportArrayDists(); // export attributes of all objects for (RooAbsArg *arg : _workspace.components()) { @@ -1939,6 +2450,12 @@ void RooJSONFactoryWSTool::exportAllObjects(JSONNode &n) _domains->writeJSON(n["domains"]); _domains.reset(); _rootnodeOutput = nullptr; + _exportedObjectNames.clear(); + _arrayComponents.clear(); + _arrayOfComponent.clear(); + _constraintTermsOfPdf.clear(); + _writtenConstraintData.clear(); + _constraintObservables.clear(); } /** @@ -2209,6 +2726,8 @@ void RooJSONFactoryWSTool::importAllNodes(const JSONNode &n) } } + warnAboutUnusedArrayComponents(_workspace); + _rootnodeInput = nullptr; _hs3VersionInput.clear(); _domains.reset(); diff --git a/roofit/hs3/test/testRooFitHS3.cxx b/roofit/hs3/test/testRooFitHS3.cxx index 560aea148a51b..750efae4457bc 100644 --- a/roofit/hs3/test/testRooFitHS3.cxx +++ b/roofit/hs3/test/testRooFitHS3.cxx @@ -9,9 +9,12 @@ #include #include #include +#include +#include #include #include #include +#include #include #include #include @@ -44,13 +47,17 @@ #include #include +#include +#include #include #include +#include #include #include #include #include +#include #include #include @@ -395,6 +402,125 @@ void importDuplicateHistoSysModel(RooWorkspace &ws, int interpCode, std::vector< ws.import(model, RooFit::Silence()); } +// A two-channel HistFactory model with a normsys shared by the channels, a staterror in both channels, a shapefactor, +// and binned data, with the constraints either written one per parameter or collected into one array distribution +// per type (`arrays`), as RooFit writes them. +std::string twoChannelHistFactoryJSON(bool arrays) +{ + const std::vector gammas{"gamma_c1_0", "gamma_c1_1", "gamma_c2_0", "gamma_c2_1", "gs_0", "gs_1"}; + std::string constraintDists; + std::string components; + std::string functions; + std::string values; + std::string domainAxes; + for (auto const &g : gammas) { + functions += R"({"name": ")" + g + R"(_poisMean", "type": "product", "factors": [")" + g + R"(", 100.0]},)"; + values += + R"({"name": "nom_)" + g + R"(", "value": 100.0, "const": true}, {"name": ")" + g + R"(", "value": 1.0},)"; + domainAxes += R"({"name": ")" + g + R"(", "min": 0.0, "max": 5.0}, {"name": "nom_)" + g + + R"(", "min": 0.0, "max": 1000.0},)"; + constraintDists += R"({"name": ")" + g + R"(_constraint", "type": "poisson_dist", "x": "nom_)" + g + + R"(", "mean": ")" + g + R"(_poisMean", "integer": false},)"; + components += R"({"x": "nom_)" + g + R"(", "mean": ")" + g + R"(_poisMean"},)"; + } + components.pop_back(); + auto constraint = [&](std::string const &individual, std::string const &array) { + return "\"" + (arrays ? array : individual) + "\""; + }; + auto gammaConstraints = [&](std::string const &prefix) { + return constraint(prefix + "_0_constraint", "pois_constraints") + ", " + + constraint(prefix + "_1_constraint", "pois_constraints"); + }; + const std::string distributions = + arrays + ? R"({"name": "gauss_constraints", "type": "array_dist", "component_type": "gaussian_dist", + "components": [{"x": "alpha_a", "mean": 0.0, "sigma": 1.0}]}, + {"name": "pois_constraints", "type": "array_dist", "component_type": "poisson_dist", + "default_parameters": {"integer": false}, "components": [)" + + components + "]}" + : R"({"name": "alpha_a_constraint", "type": "gaussian_dist", "x": "alpha_a", "mean": 0.0, "sigma": 1.0},)" + + constraintDists.substr(0, constraintDists.size() - 1); + return R"({ + "metadata": {"hs3_version": "0.3"}, + "distributions": [ + {"name": "model_c1", "type": "histfactory_dist", + "axes": [{"name": "x", "min": 0.0, "max": 2.0, "nbins": 2}], + "default_interpolation": {"type": "mult", "in": "poly6", "out": "exp"}, + "samples": [ + {"name": "sig", "data": {"contents": [10.0, 20.0], "errors": [1.0, 2.0]}, + "modifiers": [ + {"type": "normfactor", "parameters": ["mu"]}, + {"name": "a", "type": "normsys", "parameters": ["alpha_a"], "data": {"lo": 0.9, "hi": 1.1}, + "constraints": [)" + + constraint("alpha_a_constraint", "gauss_constraints") + R"(]}, + {"type": "staterror", "parameters": ["gamma_c1_0", "gamma_c1_1"], "constraints": [)" + + gammaConstraints("gamma_c1") + R"(]}]}, + {"name": "bkg", "data": {"contents": [5.0, 5.0], "errors": [1.0, 1.0]}, + "modifiers": [ + {"name": "shape", "type": "shapefactor", "parameters": ["gs_0", "gs_1"], "constraints": [)" + + gammaConstraints("gs") + R"(]}, + {"type": "staterror", "parameters": ["gamma_c1_0", "gamma_c1_1"], "constraints": [)" + + gammaConstraints("gamma_c1") + R"(]}]}]}, + {"name": "model_c2", "type": "histfactory_dist", + "axes": [{"name": "x", "min": 0.0, "max": 2.0, "nbins": 2}], + "default_interpolation": {"type": "mult", "in": "poly6", "out": "exp"}, + "samples": [ + {"name": "sig", "data": {"contents": [30.0, 10.0], "errors": [2.0, 1.0]}, + "modifiers": [ + {"type": "normfactor", "parameters": ["mu"]}, + {"name": "a", "type": "normsys", "parameters": ["alpha_a"], "data": {"lo": 0.8, "hi": 1.2}, + "constraints": [)" + + constraint("alpha_a_constraint", "gauss_constraints") + R"(]}, + {"type": "staterror", "parameters": ["gamma_c2_0", "gamma_c2_1"], "constraints": [)" + + gammaConstraints("gamma_c2") + R"(]}]}]}, + )" + + distributions + R"( + ], + "functions": [)" + + functions.substr(0, functions.size() - 1) + R"(], + "data": [ + {"name": "obs_c1", "type": "binned", "axes": [{"name": "x", "min": 0.0, "max": 2.0, "nbins": 2}], + "contents": [12.0, 19.0]}, + {"name": "obs_c2", "type": "binned", "axes": [{"name": "x", "min": 0.0, "max": 2.0, "nbins": 2}], + "contents": [28.0, 11.0]} + ], + "domains": [{"name": "default_domain", "type": "product_domain", "axes": [ + {"name": "mu", "min": 0.0, "max": 5.0}, {"name": "alpha_a", "min": -5.0, "max": 5.0}, )" + + domainAxes.substr(0, domainAxes.size() - 1) + R"(]}], + "parameter_points": [{"name": "default_values", "parameters": [ + {"name": "mu", "value": 1.0}, {"name": "alpha_a", "value": 0.0}, )" + + values.substr(0, values.size() - 1) + R"(]}], + "likelihoods": [{"name": "lh", "distributions": ["model_c1", "model_c2"], "data": ["obs_c1", "obs_c2"]}], + "analyses": [{"name": "an", "likelihood": "lh", "domains": ["default_domain"], + "parameters_of_interest": ["mu"]}] + })"; +} + +// The NLL of the two-channel model of twoChannelHistFactoryJSON() at a parameter point. +double twoChannelNLL(RooWorkspace &ws, std::map const &point) +{ + auto *mc = dynamic_cast(ws.obj("an")); + if (!mc) { + return std::numeric_limits::quiet_NaN(); + } + auto *simPdf = dynamic_cast(mc->GetPdf()); + if (!simPdf) { + return std::numeric_limits::quiet_NaN(); + } + for (auto const &[name, value] : point) { + ws.var(name)->setVal(value); + } + std::map channels; + for (auto const &item : simPdf->indexCat()) { + channels[item.first] = static_cast(ws.data("obs_" + item.first.substr(6))); + } + auto &indexCat = dynamic_cast(const_cast(simPdf->indexCat())); + RooDataHist combined{"combined", "combined", RooArgSet{*ws.var("x")}, RooFit::Index(indexCat), + RooFit::Import(channels)}; + std::unique_ptr nll{simPdf->createNLL(combined)}; + return nll->getVal(); +} + } // namespace // Test that the IO of attributes and string attributes works. @@ -514,12 +640,11 @@ TEST(RooFitHS3, ProductDomainEntriesExportExplicitBounds) EXPECT_DOUBLE_EQ((*xAxis)["min"].val_double(), -10.0); EXPECT_DOUBLE_EQ((*xAxis)["max"].val_double(), 10.0); + // Infinite bounds are expressed by omitting them, since JSON has no literal for infinity. auto const *meanAxis = RooJSONFactoryWSTool::findNamedChild((*defaultDomain)["axes"], "mean"); ASSERT_NE(meanAxis, nullptr); - ASSERT_TRUE(meanAxis->has_child("min")); - ASSERT_TRUE(meanAxis->has_child("max")); - EXPECT_TRUE((*meanAxis)["min"].is_null()); - EXPECT_TRUE((*meanAxis)["max"].is_null()); + EXPECT_FALSE(meanAxis->has_child("min")); + EXPECT_FALSE(meanAxis->has_child("max")); RooWorkspace imported; ASSERT_TRUE(RooJSONFactoryWSTool{imported}.importJSONfromString(json)); @@ -529,6 +654,21 @@ TEST(RooFitHS3, ProductDomainEntriesExportExplicitBounds) EXPECT_LT(importedMean->getMin(), 0.0); EXPECT_TRUE(std::isinf(importedMean->getMax())); EXPECT_GT(importedMean->getMax(), 0.0); + + // Files written by older RooFit versions express infinite bounds as null, which is still understood. + const std::string legacyJson = R"({ + "metadata": {"hs3_version": "0.2"}, + "domains": [{"name": "default_domain", "type": "product_domain", + "axes": [{"name": "a", "min": null, "max": 3.0}, {"name": "b", "min": 1.0}]}] + })"; + RooWorkspace legacy; + ASSERT_TRUE(RooJSONFactoryWSTool{legacy}.importJSONfromString(legacyJson)); + ASSERT_NE(legacy.var("a"), nullptr); + EXPECT_TRUE(std::isinf(legacy.var("a")->getMin())); + EXPECT_DOUBLE_EQ(legacy.var("a")->getMax(), 3.0); + ASSERT_NE(legacy.var("b"), nullptr); + EXPECT_DOUBLE_EQ(legacy.var("b")->getMin(), 1.0); + EXPECT_TRUE(std::isinf(legacy.var("b")->getMax())); } TEST(RooFitHS3, ProductDomainEntriesExportBinning) @@ -865,14 +1005,14 @@ TEST(RooFitHS3, RooGaussianConstVarSigmaExport) EXPECT_NE(domainAxes.find("\"name\":\"sigma_real\""), std::string::npos) << domainAxes; // The unbounded constant RooRealVar is still a RooRealVar, so it gets a - // domain axis with explicit null bounds that distinguishes it from a RooConstVar. + // domain axis without bounds that distinguishes it from a RooConstVar. auto tree = RooFit::Detail::JSONTree::create(json); auto const *defaultDomain = RooJSONFactoryWSTool::findNamedChild(tree->rootnode()["domains"], "default_domain"); ASSERT_NE(defaultDomain, nullptr); auto const *meanAxis = RooJSONFactoryWSTool::findNamedChild((*defaultDomain)["axes"], "mean"); ASSERT_NE(meanAxis, nullptr); - EXPECT_TRUE((*meanAxis)["min"].is_null()); - EXPECT_TRUE((*meanAxis)["max"].is_null()); + EXPECT_FALSE(meanAxis->has_child("min")); + EXPECT_FALSE(meanAxis->has_child("max")); RooWorkspace imported; ASSERT_TRUE(RooJSONFactoryWSTool{imported}.importJSONfromString(json)); @@ -2622,7 +2762,9 @@ TEST(RooFitHS3, HistFactoryConstraintKeyMigration) ASSERT_TRUE(RooJSONFactoryWSTool{ws}.importJSONfromString(jsonStr)); const std::string exported = RooJSONFactoryWSTool{ws}.exportJSONtoString(); - EXPECT_NE(exported.find("\"constraint\":\"muConstraint\""), std::string::npos) << exported; + // The constraint is collected into the array distribution of its type and top-level pdf. + EXPECT_NE(exported.find(R"("constraints":["model_channel0_gaussian_constraints"])"), std::string::npos) << exported; + EXPECT_EQ(exported.find("\"constraint\":"), std::string::npos) << exported; EXPECT_EQ(exported.find("constraint_name"), std::string::npos) << exported; EXPECT_EQ(exported.find("constraint_type"), std::string::npos) << exported; } @@ -2687,8 +2829,8 @@ TEST(RooFitHS3, HistFactoryLegacyConstraintTypeInConstraintKey) EXPECT_DOUBLE_EQ(alpha->getError(), 1.0); const std::string exported = RooJSONFactoryWSTool{ws}.exportJSONtoString(); - EXPECT_NE(exported.find("\"constraint\":\"alpha_lumiConstraint\""), std::string::npos) << exported; - EXPECT_EQ(exported.find("\"constraint\":\"Gauss\""), std::string::npos) << exported; + EXPECT_NE(exported.find(R"("constraints":["model_channel0_gaussian_constraints"])"), std::string::npos) << exported; + EXPECT_EQ(exported.find("\"constraint\":"), std::string::npos) << exported; EXPECT_EQ(exported.find("constraint_name"), std::string::npos) << exported; EXPECT_EQ(exported.find("constraint_type"), std::string::npos) << exported; } @@ -3330,9 +3472,9 @@ TEST(RooFitHS3, HistFactoryModifierWithoutConstraintIsUnconstrained) auto tree = RooFit::Detail::JSONTree::create(json1); const auto *channel = RooJSONFactoryWSTool::findNamedChild(tree->rootnode()["distributions"], "model_channel0"); ASSERT_NE(channel, nullptr) << json1; - EXPECT_FALSE(findModifier(*channel, "norm")->has_child("constraint")) << json1; - EXPECT_FALSE(findModifier(*channel, "shape")->has_child("constraint")) << json1; - EXPECT_TRUE(findModifier(*channel, "constrained")->has_child("constraint")) << json1; + EXPECT_FALSE(findModifier(*channel, "norm")->has_child("constraints")) << json1; + EXPECT_FALSE(findModifier(*channel, "shape")->has_child("constraints")) << json1; + EXPECT_TRUE(findModifier(*channel, "constrained")->has_child("constraints")) << json1; RooWorkspace ws2; ASSERT_TRUE(RooJSONFactoryWSTool{ws2}.importJSONfromString(json1)); @@ -3340,8 +3482,9 @@ TEST(RooFitHS3, HistFactoryModifierWithoutConstraintIsUnconstrained) EXPECT_EQ(RooJSONFactoryWSTool{ws2}.exportJSONtoString(), json1); } -// The constraint type of the staterror is read from the modifier, with the -// channel-level setting as fallback, and written to both places. +// The legacy constraint type of the staterror is read from the modifier, with +// the channel-level setting as fallback. On export, the constraint +// distributions are referenced by name. TEST(RooFitHS3, HistFactoryStatErrorConstraint) { RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); @@ -3386,14 +3529,23 @@ TEST(RooFitHS3, HistFactoryStatErrorConstraint) ASSERT_TRUE(RooJSONFactoryWSTool{ws}.importJSONfromString(makeJSON("Gauss", "Poisson"))); EXPECT_TRUE(hasGaussianStatConstraints(ws)); + // The constraints are exported as one array distribution that the staterror modifier references by name. const std::string json = RooJSONFactoryWSTool{ws}.exportJSONtoString(); auto tree = RooFit::Detail::JSONTree::create(json); const auto *channel = RooJSONFactoryWSTool::findNamedChild(tree->rootnode()["distributions"], "model_channel0"); ASSERT_NE(channel, nullptr) << json; - EXPECT_EQ((*channel)["staterror"]["constraint_type"].val(), "Gauss") << json; + EXPECT_FALSE(channel->has_child("staterror")) << json; const auto *modifier = findModifier(*channel, "staterror"); ASSERT_NE(modifier, nullptr) << json; - EXPECT_EQ((*modifier)["constraint"].val(), "Gauss") << json; + EXPECT_FALSE(modifier->has_child("constraint")) << json; + ASSERT_EQ((*modifier)["constraints"].num_children(), 2u) << json; + for (auto const &constraintName : (*modifier)["constraints"].children()) { + const auto *constraint = + RooJSONFactoryWSTool::findNamedChild(tree->rootnode()["distributions"], constraintName.val()); + ASSERT_NE(constraint, nullptr) << json; + EXPECT_EQ((*constraint)["type"].val(), "array_dist") << json; + EXPECT_EQ((*constraint)["component_type"].val(), "gaussian_dist") << json; + } RooWorkspace ws2; ASSERT_TRUE(RooJSONFactoryWSTool{ws2}.importJSONfromString(json)); @@ -3443,12 +3595,9 @@ TEST(RooFitHS3, HistFactoryMissingInterpolation) // RooFit always writes a "default_interpolation" if a modifier relies on it. const std::string exported = RooJSONFactoryWSTool{ws}.exportJSONtoString(); EXPECT_NE(exported.find("default_interpolation"), std::string::npos) << exported; + EXPECT_NE(exported.find(R"("hs3_version":"0.3")"), std::string::npos) << exported; RooWorkspace ws2; - std::string exportedCurrent = exported; - const std::string exportedVersion = R"("hs3_version":"0.2")"; - ASSERT_NE(exportedCurrent.find(exportedVersion), std::string::npos) << exported; - exportedCurrent.replace(exportedCurrent.find(exportedVersion), exportedVersion.size(), R"("hs3_version":"0.3")"); - EXPECT_TRUE(RooJSONFactoryWSTool{ws2}.importJSONfromString(exportedCurrent)); + EXPECT_TRUE(RooJSONFactoryWSTool{ws2}.importJSONfromString(exported)); } // Malformed matrices, error arrays and bin edges are rejected instead of being read out of bounds. @@ -3516,6 +3665,1022 @@ TEST(RooFitHS3, UnbinnedAxisRangesOnImport) } } +// The analysis keys "parameters_of_interest" and "domains" are written. The +// singular keys of HS3 v0.2.9 are still read, with a deprecation warning. +TEST(RooFitHS3, AnalysisLegacySingularKeys) +{ + RooWorkspace ws1; + ws1.factory("Gaussian::pdf(x[0, -5, 5], mu[0, -5, 5], sigma[1, 0.1, 5])"); + RooFit::ModelConfig mc1{"ModelConfig", &ws1}; + mc1.SetPdf("pdf"); + mc1.SetParametersOfInterest("mu"); + mc1.SetNuisanceParameters("sigma"); + mc1.SetObservables("x"); + ws1.import(mc1); + + const std::string json = RooJSONFactoryWSTool{ws1}.exportJSONtoString(); + EXPECT_NE(json.find(R"("parameters_of_interest":["mu"])"), std::string::npos) << json; + EXPECT_NE(json.find(R"("domains":["ModelConfig_nuisance_parameters","ModelConfig_parameters_of_interest"])"), + std::string::npos) + << json; + + // Turn the analysis into the singular HS3 v0.2.9 form, with only the nuisance parameter domain. + std::string legacyJson = json; + auto replace = [&](std::string const &from, std::string const &to) { + const auto pos = legacyJson.find(from); + ASSERT_NE(pos, std::string::npos) << legacyJson; + legacyJson.replace(pos, from.size(), to); + }; + replace(R"("parameters_of_interest":["mu"])", R"("parameter_of_interest":["mu"])"); + replace(R"("domains":["ModelConfig_nuisance_parameters","ModelConfig_parameters_of_interest"])", + R"("domain":"ModelConfig_nuisance_parameters")"); + + RooWorkspace ws2; + const std::string warnings = captureMessages(RooFit::WARNING, RooFit::ObjectHandling, [&] { + EXPECT_TRUE(RooJSONFactoryWSTool{ws2}.importJSONfromString(legacyJson)); + }); + EXPECT_NE(warnings.find(R"(the key "parameter_of_interest" of analysis)"), std::string::npos) << warnings; + EXPECT_NE(warnings.find(R"(the key "domain" of analysis)"), std::string::npos) << warnings; + auto *mc2 = dynamic_cast(ws2.obj("ModelConfig")); + ASSERT_NE(mc2, nullptr); + EXPECT_NE(mc2->GetParametersOfInterest()->find("mu"), nullptr); + EXPECT_NE(mc2->GetNuisanceParameters()->find("sigma"), nullptr); +} + +// "crystalball_dist" is the single-sided Crystal Ball, "crystalball_doublesided_dist" the double-sided one. Before +// HS3 0.3, "crystalball_dist" with the "_L"/"_R" fields was the double-sided one, which is still read. +TEST(RooFitHS3, CrystalBallTypes) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + RooRealVar x{"x", "x", 0.0, 10.0}; + x.setBins(20); + RooRealVar mean{"mean", "mean", 5.0, 0.0, 10.0}; + RooRealVar sigma{"sigma", "sigma", 1.0, 0.1, 5.0}; + RooRealVar sigmaR{"sigmaR", "sigmaR", 1.5, 0.1, 5.0}; + RooRealVar alphaL{"alphaL", "alphaL", 1.2, 0.1, 5.0}; + RooRealVar nL{"nL", "nL", 3.0, 0.5, 10.0}; + RooRealVar alphaR{"alphaR", "alphaR", 2.0, 0.1, 5.0}; + RooRealVar nR{"nR", "nR", 5.0, 0.5, 10.0}; + RooRealVar alphaNeg{"alphaNeg", "alphaNeg", -1.2, -5.0, -0.1}; + + auto exportToString = [](RooAbsArg const &arg) { + RooWorkspace ws; + ws.import(arg, RooFit::Silence()); + return RooJSONFactoryWSTool{ws}.exportJSONtoString(); + }; + + // RooCBShape is the single-sided crystalball_dist. + RooCBShape cbShape{"cbShape", "", x, mean, sigma, alphaL, nL}; + EXPECT_EQ(validate(cbShape), 0); + EXPECT_NE(exportToString(cbShape).find(R"("type":"crystalball_dist")"), std::string::npos); + + // A single-sided RooCrystalBall is a crystalball_dist too, for both signs of alpha. It is read back as a + // RooCBShape, which is the same function. + for (RooRealVar *alpha : {&alphaL, &alphaNeg}) { + RooCrystalBall cbSingle{"cbSingle", "", x, mean, sigma, *alpha, nL, /*doubleSided=*/false}; + const std::string json = exportToString(cbSingle); + EXPECT_NE(json.find(R"("type":"crystalball_dist")"), std::string::npos) << json; + EXPECT_NE(json.find(std::string(R"("alpha":")") + alpha->GetName() + "\""), std::string::npos) << json; + EXPECT_EQ(json.find("alpha_L"), std::string::npos) << json; + + // The class changes on the way back, so validate() can't be used: the parameter points are written in a + // different order. Compare the values instead, and check that the JSON is stable from the second export on. + RooWorkspace ws2; + RooJSONFactoryWSTool{ws2}.importJSONfromString(json); + auto *cbShapeBack = dynamic_cast(ws2.pdf("cbSingle")); + ASSERT_NE(cbShapeBack, nullptr); + RooRealVar &x2 = *ws2.var("x"); + RooArgSet nset1{x}; + RooArgSet nset2{x2}; + for (int i = 0; i < x.numBins(); ++i) { + x.setBin(i); + x2.setBin(i); + EXPECT_NEAR(cbShapeBack->getVal(nset2), cbSingle.getVal(nset1), 1e-10) << alpha->GetName() << " bin " << i; + } + const std::string json2 = RooJSONFactoryWSTool{ws2}.exportJSONtoString(); + RooWorkspace ws3; + RooJSONFactoryWSTool{ws3}.importJSONfromString(json2); + EXPECT_EQ(RooJSONFactoryWSTool{ws3}.exportJSONtoString(), json2); + } + + // Double-sided RooCrystalBall shapes are crystalball_doublesided_dist, with the tails always written for both + // sides, and the width written as "sigma" if both sides share it. + RooCrystalBall cbSymmetric{"cbSymmetric", "", x, mean, sigma, alphaL, nL, /*doubleSided=*/true}; + RooCrystalBall cbAsymmetric{"cbAsymmetric", "", x, mean, sigma, sigmaR, alphaL, nL, alphaR, nR}; + EXPECT_EQ(validate(cbSymmetric), 0); + EXPECT_EQ(validate(cbAsymmetric), 0); + { + const std::string json = exportToString(cbSymmetric); + EXPECT_NE(json.find(R"("type":"crystalball_doublesided_dist")"), std::string::npos) << json; + EXPECT_NE(json.find(R"("sigma":"sigma")"), std::string::npos) << json; + EXPECT_NE(json.find(R"("alpha_R":"alphaL")"), std::string::npos) << json; + EXPECT_EQ(json.find(R"("alpha":)"), std::string::npos) << json; + } + { + const std::string json = exportToString(cbAsymmetric); + EXPECT_NE(json.find(R"("type":"crystalball_doublesided_dist")"), std::string::npos) << json; + EXPECT_NE(json.find(R"("sigma_R":"sigmaR")"), std::string::npos) << json; + EXPECT_NE(json.find(R"("alpha_R":"alphaR")"), std::string::npos) << json; + } + + auto document = [](std::string const &distribution) { + return R"({ + "metadata": {"hs3_version": "0.3"}, + "domains": [{"name": "default_domain", "type": "product_domain", "axes": [{"name": "m", "min": 0, "max": 10}]}], + "distributions": [)" + + distribution + R"(] + })"; + }; + auto importPdf = [&](std::string const &distribution, std::string *warnings = nullptr) { + auto ws = std::make_unique(); + const std::string messages = captureMessages(RooFit::WARNING, RooFit::ObjectHandling, [&] { + RooJSONFactoryWSTool{*ws}.importJSONfromString(document(distribution)); + }); + if (warnings) { + *warnings = messages; + } + return ws; + }; + + // Import of the two standard types. + auto single = importPdf(R"({"name": "cb", "type": "crystalball_dist", "m": "m", "m0": 5.0, "sigma": 1.0, + "alpha": -1.2, "n": 3.0})"); + EXPECT_NE(dynamic_cast(single->pdf("cb")), nullptr); + auto doubleSided = importPdf(R"({"name": "cb", "type": "crystalball_doublesided_dist", "m": "m", "m0": 5.0, + "sigma": 1.0, "alpha_L": 1.2, "n_L": 3.0, "alpha_R": 2.0, "n_R": 5.0})"); + EXPECT_NE(dynamic_cast(doubleSided->pdf("cb")), nullptr); + + // The double-sided form of crystalball_dist from before HS3 0.3 is still read, with a deprecation warning. This + // includes the symmetric tails given as "alpha"/"n" with asymmetric widths. + std::string warnings; + auto legacy = importPdf(R"({"name": "cb", "type": "crystalball_dist", "m": "m", "m0": 5.0, "sigma": 1.0, + "alpha_L": 1.2, "n_L": 3.0, "alpha_R": 2.0, "n_R": 5.0})", + &warnings); + auto *legacyCb = dynamic_cast(legacy->pdf("cb")); + ASSERT_NE(legacyCb, nullptr); + ASSERT_NE(legacyCb->alphaR(), nullptr); + EXPECT_DOUBLE_EQ(legacyCb->alphaR()->getVal(), 2.0); + EXPECT_NE(warnings.find("please use \"crystalball_doublesided_dist\""), std::string::npos) << warnings; + auto legacySymmetric = importPdf(R"({"name": "cb", "type": "crystalball_dist", "m": "m", "m0": 5.0, + "sigma_L": 1.0, "sigma_R": 1.5, "alpha": 1.2, "n": 3.0})"); + auto *legacySymmetricCb = dynamic_cast(legacySymmetric->pdf("cb")); + ASSERT_NE(legacySymmetricCb, nullptr); + ASSERT_NE(legacySymmetricCb->alphaR(), nullptr); + EXPECT_DOUBLE_EQ(legacySymmetricCb->alphaR()->getVal(), 1.2); + + // Invalid field combinations. + expectImportThrowsWithError(document(R"({"name": "cb", "type": "crystalball_doublesided_dist", "m": "m", + "m0": 5.0, "sigma": 1.0, "alpha": 1.2, "n": 3.0})"), + "must not have \"alpha\" or \"n\""); + expectImportThrowsWithError(document(R"({"name": "cb", "type": "crystalball_doublesided_dist", "m": "m", + "m0": 5.0, "sigma": 1.0, "alpha_L": 1.2, "n_L": 3.0, "alpha_R": 2.0})"), + "is missing \"n_R\""); + expectImportThrowsWithError(document(R"({"name": "cb", "type": "crystalball_doublesided_dist", "m": "m", + "m0": 5.0, "sigma": 1.0, "sigma_L": 1.0, "sigma_R": 1.0, "alpha_L": 1.2, "n_L": 3.0, "alpha_R": 2.0, + "n_R": 5.0})"), + "must not combine \"sigma\""); + expectImportThrowsWithError(document(R"({"name": "cb", "type": "crystalball_doublesided_dist", "m": "m", + "m0": 5.0, "sigma_L": 1.0, "alpha_L": 1.2, "n_L": 3.0, "alpha_R": 2.0, "n_R": 5.0})"), + "needs both \"sigma_L\" and \"sigma_R\""); + expectImportThrowsWithError(document(R"({"name": "cb", "type": "crystalball_dist", "m": "m", "m0": 5.0, + "sigma": 1.0, "alpha": 1.2, "n": 3.0, "alpha_R": 2.0, "n_R": 5.0})"), + "must not combine"); +} + +// HistFactory modifiers reference their parameters with "parameters" and their constraint distributions by name with +// "constraints", where null leaves a parameter unconstrained. +TEST(RooFitHS3, HistFactoryConstraintsByName) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + const std::string json = R"({ + "metadata": {"hs3_version": "0.3"}, + "distributions": [ + { + "name": "model_channel0", + "type": "histfactory_dist", + "axes": [{"name": "x", "min": 0.0, "max": 2.0, "nbins": 2}], + "default_interpolation": {"type": "mult", "in": "poly6", "out": "exp"}, + "samples": [ + { + "name": "sig", + "data": {"contents": [10.0, 20.0], "errors": [1.0, 2.0]}, + "modifiers": [ + {"type": "normfactor", "parameters": ["mu"]}, + {"type": "normsys", "parameters": ["alpha_a"], "constraints": ["alpha_a_constraint"], + "data": {"lo": 0.9, "hi": 1.1}}, + {"type": "normsys", "parameters": ["alpha_b"], "constraints": [null], + "data": {"lo": 0.8, "hi": 1.2}}, + {"type": "staterror", "parameters": ["gamma_0", "gamma_1"], + "constraints": ["gamma_0_constraint", null]} + ] + } + ] + }, + {"name": "alpha_a_constraint", "type": "gaussian_dist", "x": "alpha_a", "mean": 0.0, "sigma": 1.0}, + {"name": "gamma_0_constraint", "type": "gaussian_dist", "x": "gamma_0", "mean": 1.0, "sigma": 0.1} + ], + "domains": [{"name": "default_domain", "type": "product_domain", "axes": [ + {"name": "alpha_a", "min": -5, "max": 5}, {"name": "alpha_b", "min": -5, "max": 5}, + {"name": "gamma_0", "min": 0, "max": 5}, {"name": "gamma_1", "min": 0, "max": 5}]}] + })"; + + RooWorkspace ws; + ASSERT_TRUE(RooJSONFactoryWSTool{ws}.importJSONfromString(json)); + auto *model = dynamic_cast(ws.pdf("model_channel0")); + ASSERT_NE(model, nullptr); + RooArgSet constraintNames; + for (RooAbsArg *pdf : model->pdfList()) { + if (!dynamic_cast(pdf)) { + constraintNames.add(*pdf); + } + } + // Exactly the two named constraints, and no constraints created by RooFit. + EXPECT_EQ(constraintNames.size(), 2u); + EXPECT_NE(constraintNames.find("alpha_a_constraint"), nullptr); + EXPECT_NE(constraintNames.find("gamma_0_constraint"), nullptr); + EXPECT_NE(ws.var("mu"), nullptr); + + // The export collects the constraints into the array distribution of their type (both are Gaussians), the + // unconstrained parameter keeps its null, and the result is stable under a round trip. + const std::string exported = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + EXPECT_NE(exported.find(R"("constraints":["model_channel0_gaussian_constraints"])"), std::string::npos) << exported; + EXPECT_NE(exported.find(R"("constraints":["model_channel0_gaussian_constraints",null])"), std::string::npos) + << exported; + EXPECT_EQ(exported.find(R"("parameter":)"), std::string::npos) << exported; + RooWorkspace ws2; + ASSERT_TRUE(RooJSONFactoryWSTool{ws2}.importJSONfromString(exported)); + EXPECT_EQ(RooJSONFactoryWSTool{ws2}.exportJSONtoString(), exported); + + // In HS3 >= 0.3, a staterror without any constraint information is unconstrained. + std::string unconstrained = json; + const std::string statConstraints = R"(, + "constraints": ["gamma_0_constraint", null])"; + ASSERT_NE(unconstrained.find(statConstraints), std::string::npos); + unconstrained.replace(unconstrained.find(statConstraints), statConstraints.size(), ""); + RooWorkspace ws3; + ASSERT_TRUE(RooJSONFactoryWSTool{ws3}.importJSONfromString(unconstrained)); + auto *model3 = dynamic_cast(ws3.pdf("model_channel0")); + ASSERT_NE(model3, nullptr); + EXPECT_EQ(model3->pdfList().size(), 2u); // the sum and alpha_a_constraint +} + +// An array_dist is imported as one distribution per component, named after its observable, held together by a +// RooProdPdf with the attribute hs3_array_dist that carries the name of the array. The component types are built by +// their own importers, with "default_parameters" filled in where a component omits a parameter. +TEST(RooFitHS3, ArrayDistImport) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + const std::string json = R"({ + "metadata": {"hs3_version": "0.3"}, + "distributions": [ + {"name": "pois", "type": "array_dist", "component_type": "poisson_dist", + "default_parameters": {"mean": "lambda_default"}, + "components": [{"x": "n_1", "mean": "lambda_1"}, {"x": "n_2"}, {"x": "n_3", "mean": 7.5, "integer": true}]}, + {"name": "gauss", "type": "array_dist", "component_type": "gaussian_dist", + "default_parameters": {"sigma": 1.0}, + "components": [{"x": "nom_a", "mean": "a", "sigma": "sigma_a"}, {"x": "nom_b", "mean": "b"}, + {"x": "nom_c", "mean": 0.5, "sigma": 2.0}]}, + {"name": "logn", "type": "array_dist", "component_type": "lognormal_dist", + "components": [{"x": "y_1", "mu": "mu_1", "sigma": 1.5}, {"x": "y_2", "mu": 2.0, "sigma": "sigma_y"}]} + ], + "functions": [ + {"name": "lambda_1", "type": "product", "factors": [20.0, "np_1"]}, + {"name": "lambda_default", "type": "product", "factors": [10.0, "np_2"]} + ], + "domains": [{"name": "default_domain", "type": "product_domain", "axes": [ + {"name": "n_1", "min": 0.0, "max": 100.0}, {"name": "n_2", "min": 0.0, "max": 100.0}, + {"name": "n_3", "min": 0.0, "max": 100.0}, + {"name": "np_1", "min": 0.0, "max": 5.0}, {"name": "np_2", "min": 0.0, "max": 5.0}, + {"name": "nom_a", "min": -5.0, "max": 5.0}, {"name": "nom_b", "min": -5.0, "max": 5.0}, + {"name": "nom_c", "min": -5.0, "max": 5.0}, {"name": "a", "min": -5.0, "max": 5.0}, + {"name": "b", "min": -5.0, "max": 5.0}, {"name": "sigma_a", "min": 0.1, "max": 10.0}, + {"name": "y_1", "min": 0.1, "max": 10.0}, {"name": "y_2", "min": 0.1, "max": 10.0}, + {"name": "mu_1", "min": 0.1, "max": 10.0}, {"name": "sigma_y", "min": 1.0, "max": 10.0}]}], + "parameter_points": [{"name": "default_values", "parameters": [ + {"name": "n_1", "value": 20.0}, {"name": "n_2", "value": 10.0}, {"name": "n_3", "value": 8.0}, + {"name": "np_1", "value": 1.1}, {"name": "np_2", "value": 0.9}, + {"name": "nom_a", "value": 0.3}, {"name": "nom_b", "value": -0.2}, {"name": "nom_c", "value": 1.0}, + {"name": "a", "value": 0.1}, {"name": "b", "value": 0.4}, {"name": "sigma_a", "value": 0.5}, + {"name": "y_1", "value": 2.0}, {"name": "y_2", "value": 3.0}, {"name": "mu_1", "value": 1.5}, + {"name": "sigma_y", "value": 1.2}]}] + })"; + + RooWorkspace ws; + ASSERT_TRUE(RooJSONFactoryWSTool{ws}.importJSONfromString(json)); + + auto checkArray = [&](const char *name, std::vector const &components) { + auto *array = dynamic_cast(ws.pdf(name)); + ASSERT_NE(array, nullptr) << name; + EXPECT_TRUE(array->getAttribute("hs3_array_dist")) << name; + ASSERT_EQ(array->pdfList().size(), components.size()) << name; + double product = 1.0; + for (std::size_t i = 0; i < components.size(); ++i) { + EXPECT_STREQ(array->pdfList()[i].GetName(), components[i].c_str()) << name; + product *= static_cast(array->pdfList()[i]).getVal(); + } + // The array is the product of its components. + EXPECT_DOUBLE_EQ(array->getVal(), product) << name; + }; + checkArray("pois", {"n_1_constraint", "n_2_constraint", "n_3_constraint"}); + checkArray("gauss", {"nom_a_constraint", "nom_b_constraint", "nom_c_constraint"}); + checkArray("logn", {"y_1_constraint", "y_2_constraint"}); + + // Poisson: own parameter, default parameter, number. Without "integer", the continued form is used. + auto *p1 = dynamic_cast(ws.pdf("n_1_constraint")); + auto *p2 = dynamic_cast(ws.pdf("n_2_constraint")); + auto *p3 = dynamic_cast(ws.pdf("n_3_constraint")); + ASSERT_NE(p1, nullptr); + ASSERT_NE(p2, nullptr); + ASSERT_NE(p3, nullptr); + EXPECT_STREQ(p1->getX().GetName(), "n_1"); + EXPECT_STREQ(p1->getMean().GetName(), "lambda_1"); + EXPECT_STREQ(p2->getMean().GetName(), "lambda_default"); + EXPECT_DOUBLE_EQ(p3->getMean().getVal(), 7.5); + EXPECT_TRUE(p1->getNoRounding()); + EXPECT_TRUE(p2->getNoRounding()); + EXPECT_FALSE(p3->getNoRounding()); + EXPECT_DOUBLE_EQ(p1->getVal(), TMath::Poisson(20.0, 22.0)); + EXPECT_DOUBLE_EQ(p2->getVal(), TMath::Poisson(10.0, 9.0)); + + // Gaussian: parameters by name, default sigma, numbers. + auto *ga = dynamic_cast(ws.pdf("nom_a_constraint")); + auto *gb = dynamic_cast(ws.pdf("nom_b_constraint")); + auto *gc = dynamic_cast(ws.pdf("nom_c_constraint")); + ASSERT_NE(ga, nullptr); + ASSERT_NE(gb, nullptr); + ASSERT_NE(gc, nullptr); + EXPECT_STREQ(ga->getX().GetName(), "nom_a"); + EXPECT_STREQ(ga->getMean().GetName(), "a"); + EXPECT_STREQ(ga->getSigma().GetName(), "sigma_a"); + EXPECT_STREQ(gb->getMean().GetName(), "b"); + EXPECT_DOUBLE_EQ(gb->getSigma().getVal(), 1.0); + EXPECT_DOUBLE_EQ(gc->getMean().getVal(), 0.5); + EXPECT_DOUBLE_EQ(gc->getSigma().getVal(), 2.0); + EXPECT_DOUBLE_EQ(ga->getVal(), std::exp(-0.5 * std::pow((0.3 - 0.1) / 0.5, 2))); + EXPECT_DOUBLE_EQ(gb->getVal(), std::exp(-0.5 * std::pow((-0.2 - 0.4) / 1.0, 2))); + EXPECT_DOUBLE_EQ(gc->getVal(), std::exp(-0.5 * std::pow((1.0 - 0.5) / 2.0, 2))); + + // LogNormal in the standard parametrization. + auto *l1 = dynamic_cast(ws.pdf("y_1_constraint")); + auto *l2 = dynamic_cast(ws.pdf("y_2_constraint")); + ASSERT_NE(l1, nullptr); + ASSERT_NE(l2, nullptr); + EXPECT_TRUE(l1->useStandardParametrization()); + EXPECT_STREQ(l1->getMedian().GetName(), "mu_1"); + EXPECT_DOUBLE_EQ(l1->getShapeK().getVal(), 1.5); + EXPECT_DOUBLE_EQ(l2->getMedian().getVal(), 2.0); + EXPECT_STREQ(l2->getShapeK().GetName(), "sigma_y"); + RooLognormal ref{"ref", "ref", *ws.var("y_1"), *ws.var("mu_1"), RooFit::RooConst(1.5), true}; + EXPECT_DOUBLE_EQ(l1->getVal(), ref.getVal()); +} + +// Invalid array distributions are rejected with explicit errors. +TEST(RooFitHS3, ArrayDistImportErrors) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + auto document = [](std::string const &array, std::string const &functions = "") { + return R"({"metadata": {"hs3_version": "0.3"}, "distributions": [)" + array + R"(], + "functions": [)" + + functions + R"(], + "domains": [{"name": "default_domain", "type": "product_domain", "axes": [ + {"name": "n_1", "min": 0.0, "max": 100.0}, {"name": "n_2", "min": 0.0, "max": 100.0}, + {"name": "mu_1", "min": 0.0, "max": 100.0}, {"name": "mu_2", "min": 0.0, "max": 100.0}]}]})"; + }; + auto array = [](std::string const &body) { return R"({"name": "arr", "type": "array_dist", )" + body + "}"; }; + + expectImportThrowsWithError(document(array(R"("components": [{"x": "n_1", "mean": "mu_1"}])")), + "missing \"component_type\""); + expectImportThrowsWithError( + document(array(R"("component_type": "exponential_dist", "components": [{"x": "n_1", "c": 1.0}])")), + "the component_type 'exponential_dist' is not supported"); + expectImportThrowsWithError(document(array(R"("component_type": "poisson_dist", "components": [])")), + "\"components\" must be a non-empty list"); + expectImportThrowsWithError( + document(array(R"("component_type": "poisson_dist", "default_parameters": {"x": "n_1", "mean": "mu_1"}, + "components": [{"x": "n_1"}])")), + "\"x\" must be given per component, not in \"default_parameters\""); + expectImportThrowsWithError(document(array(R"("component_type": "poisson_dist", + "components": [{"x": "n_1", "mean": "mu_1"}, {"x": "n_1", "mean": "mu_2"}])")), + "the observable 'n_1' is used by two components"); + expectImportThrowsWithError( + document(array(R"("component_type": "poisson_dist", "components": [{"x": 3.0, "mean": "mu_1"}])")), + "\"x\" must be the name of the observable"); + expectImportThrowsWithError(document(array(R"("component_type": "poisson_dist", "components": [{"mean": "mu_1"}])")), + "component 0 has no \"x\""); + // A parameter given neither in the component nor in the defaults: the error of the component type's importer. + expectImportThrowsWithError(document(array(R"("component_type": "poisson_dist", "components": [{"x": "n_1"}])")), + "no \"mean\" given in \"n_1_constraint\""); + expectImportThrowsWithError(document(array(R"("component_type": "poisson_dist", + "components": [{"x": "n_1", "mean": "n_2"}, {"x": "n_2", "mean": "mu_2"}])")), + "component 0 depends on the observable 'n_2' of component 1"); + expectImportThrowsWithError( + document(array(R"("component_type": "poisson_dist", "components": [{"x": "n_1", "mean": "f"}])"), + R"({"name": "f", "type": "product", "factors": ["n_1", 2.0]})"), + "the parameters of component 0 depend on its own observable 'n_1'"); +} + +// A RooProdPdf with the attribute hs3_array_dist is written as array_dist: one component per factor with the fields +// of the standalone distribution, values shared by all components in "default_parameters". The written document +// has no standalone distribution for the components, and it is stable under a round trip. +TEST(RooFitHS3, ArrayDistExport) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + const std::string json = R"({ + "metadata": {"hs3_version": "0.3"}, + "distributions": [ + {"name": "pois", "type": "array_dist", "component_type": "poisson_dist", + "default_parameters": {"mean": "lambda_default"}, + "components": [{"x": "n_1", "mean": "lambda_1"}, {"x": "n_2"}, {"x": "n_3", "mean": 7.5, "integer": true}]}, + {"name": "gauss", "type": "array_dist", "component_type": "gaussian_dist", + "default_parameters": {"sigma": 1.0}, + "components": [{"x": "nom_a", "mean": "a"}, {"x": "nom_b", "mean": "b"}, {"x": "nom_c", "mean": 0.5}]} + ], + "functions": [ + {"name": "lambda_1", "type": "product", "factors": [20.0, "np_1"]}, + {"name": "lambda_default", "type": "product", "factors": [10.0, "np_2"]} + ], + "domains": [{"name": "default_domain", "type": "product_domain", "axes": [ + {"name": "n_1", "min": 0.0, "max": 100.0}, {"name": "n_2", "min": 0.0, "max": 100.0}, + {"name": "n_3", "min": 0.0, "max": 100.0}, + {"name": "np_1", "min": 0.0, "max": 5.0}, {"name": "np_2", "min": 0.0, "max": 5.0}, + {"name": "nom_a", "min": -5.0, "max": 5.0}, {"name": "nom_b", "min": -5.0, "max": 5.0}, + {"name": "nom_c", "min": -5.0, "max": 5.0}, {"name": "a", "min": -5.0, "max": 5.0}, + {"name": "b", "min": -5.0, "max": 5.0}]}] + })"; + + RooWorkspace ws; + ASSERT_TRUE(RooJSONFactoryWSTool{ws}.importJSONfromString(json)); + const std::string exported = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + + auto tree = RooFit::Detail::JSONTree::create(exported); + auto const &distributions = tree->rootnode()["distributions"]; + // The components are not written on their own. + EXPECT_EQ(distributions.num_children(), 2u) << exported; + const auto *pois = RooJSONFactoryWSTool::findNamedChild(distributions, "pois"); + ASSERT_NE(pois, nullptr) << exported; + EXPECT_EQ((*pois)["type"].val(), "array_dist"); + EXPECT_EQ((*pois)["component_type"].val(), "poisson_dist"); + // The Poisson components have different "integer" flags, so nothing is shared. + EXPECT_FALSE(pois->has_child("default_parameters")) << exported; + ASSERT_EQ((*pois)["components"].num_children(), 3u); + EXPECT_EQ((*pois)["components"].child(0)["x"].val(), "n_1"); + EXPECT_EQ((*pois)["components"].child(0)["mean"].val(), "lambda_1"); + EXPECT_FALSE((*pois)["components"].child(0)["integer"].val_bool()); + EXPECT_EQ((*pois)["components"].child(1)["x"].val(), "n_2"); + EXPECT_EQ((*pois)["components"].child(1)["mean"].val(), "lambda_default"); + EXPECT_EQ((*pois)["components"].child(2)["x"].val(), "n_3"); + EXPECT_DOUBLE_EQ((*pois)["components"].child(2)["mean"].val_double(), 7.5); + EXPECT_TRUE((*pois)["components"].child(2)["integer"].val_bool()); + + const auto *gauss = RooJSONFactoryWSTool::findNamedChild(distributions, "gauss"); + ASSERT_NE(gauss, nullptr) << exported; + EXPECT_EQ((*gauss)["component_type"].val(), "gaussian_dist"); + ASSERT_TRUE(gauss->has_child("default_parameters")) << exported; + EXPECT_DOUBLE_EQ((*gauss)["default_parameters"]["sigma"].val_double(), 1.0); + ASSERT_EQ((*gauss)["components"].num_children(), 3u); + for (auto const &component : (*gauss)["components"].children()) { + EXPECT_FALSE(component.has_child("sigma")) << exported; + } + EXPECT_EQ((*gauss)["components"].child(0)["mean"].val(), "a"); + EXPECT_DOUBLE_EQ((*gauss)["components"].child(2)["mean"].val_double(), 0.5); + + // The tag is not written as an attribute; it is restored from the type on import. + EXPECT_EQ(exported.find("hs3_array_dist"), std::string::npos) << exported; + + RooWorkspace ws2; + ASSERT_TRUE(RooJSONFactoryWSTool{ws2}.importJSONfromString(exported)); + EXPECT_TRUE(ws2.pdf("gauss")->getAttribute("hs3_array_dist")); + EXPECT_EQ(RooJSONFactoryWSTool{ws2}.exportJSONtoString(), exported); +} + +// An array distribution built in RooFit (a tagged RooProdPdf), the opt-out that writes it as product_dist, and the +// products that cannot be written as array_dist. +TEST(RooFitHS3, ArrayDistExportOfTaggedProduct) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + RooRealVar nomA{"nom_a", "nom_a", 0.0, -5.0, 5.0}; + RooRealVar nomB{"nom_b", "nom_b", 0.0, -5.0, 5.0}; + RooRealVar a{"a", "a", 0.0, -5.0, 5.0}; + RooRealVar b{"b", "b", 0.0, -5.0, 5.0}; + RooGaussian gaussA{"gauss_a", "gauss_a", nomA, a, RooFit::RooConst(1.0)}; + RooGaussian gaussB{"gauss_b", "gauss_b", nomB, b, RooFit::RooConst(1.0)}; + RooProdPdf array{"constraints", "constraints", {gaussA, gaussB}}; + array.setAttribute("hs3_array_dist"); + + RooWorkspace ws; + ws.import(array, RooFit::Silence()); + + const std::string exported = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + auto tree = RooFit::Detail::JSONTree::create(exported); + auto const &distributions = tree->rootnode()["distributions"]; + ASSERT_EQ(distributions.num_children(), 1u) << exported; + const auto *node = RooJSONFactoryWSTool::findNamedChild(distributions, "constraints"); + ASSERT_NE(node, nullptr) << exported; + EXPECT_EQ((*node)["type"].val(), "array_dist"); + EXPECT_EQ((*node)["component_type"].val(), "gaussian_dist"); + EXPECT_DOUBLE_EQ((*node)["default_parameters"]["sigma"].val_double(), 1.0); + ASSERT_EQ((*node)["components"].num_children(), 2u); + EXPECT_EQ((*node)["components"].child(0)["x"].val(), "nom_a"); + EXPECT_EQ((*node)["components"].child(0)["mean"].val(), "a"); + EXPECT_EQ((*node)["components"].child(1)["x"].val(), "nom_b"); + + // Round trip: the components get the names of their observables, the product is tagged again. + RooWorkspace ws2; + ASSERT_TRUE(RooJSONFactoryWSTool{ws2}.importJSONfromString(exported)); + auto *imported = dynamic_cast(ws2.pdf("constraints")); + ASSERT_NE(imported, nullptr); + EXPECT_TRUE(imported->getAttribute("hs3_array_dist")); + EXPECT_NE(dynamic_cast(ws2.pdf("nom_a_constraint")), nullptr); + EXPECT_EQ(RooJSONFactoryWSTool{ws2}.exportJSONtoString(), exported); + + // Opt-out: no array_dist is written, the tagged product is an ordinary product_dist. + RooJSONFactoryWSTool::config().exportArrayDist = false; + const std::string plain = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + RooJSONFactoryWSTool::config().exportArrayDist = true; + auto plainTree = RooFit::Detail::JSONTree::create(plain); + auto const &plainDistributions = plainTree->rootnode()["distributions"]; + EXPECT_EQ(plainDistributions.num_children(), 3u) << plain; + const auto *product = RooJSONFactoryWSTool::findNamedChild(plainDistributions, "constraints"); + ASSERT_NE(product, nullptr) << plain; + EXPECT_EQ((*product)["type"].val(), "product_dist"); + const auto *component = RooJSONFactoryWSTool::findNamedChild(plainDistributions, "gauss_a"); + ASSERT_NE(component, nullptr) << plain; + EXPECT_EQ((*component)["type"].val(), "gaussian_dist"); + EXPECT_EQ(plain.find("array_dist"), std::string::npos) << plain; + + // Products that cannot be written as array_dist: a factor of another class, two factors over the same + // observable, and a tag on something that is not a product. + { + RooRealVar x{"x", "x", 1.0, 0.0, 10.0}; + RooRealVar c{"c", "c", -1.0, -10.0, 0.0}; + RooExponential expo{"expo", "expo", x, c}; + RooProdPdf mixed{"mixed", "mixed", {gaussA, expo}}; + mixed.setAttribute("hs3_array_dist"); + RooWorkspace wsMixed; + wsMixed.import(mixed, RooFit::Silence()); + expectExportThrowsWithError(wsMixed, "'expo' cannot be a component"); + } + { + RooGaussian gaussA2{"gauss_a2", "gauss_a2", nomA, b, RooFit::RooConst(2.0)}; + RooProdPdf sameX{"same_x", "same_x", {gaussA, gaussA2}}; + sameX.setAttribute("hs3_array_dist"); + RooWorkspace wsSameX; + wsSameX.import(sameX, RooFit::Silence()); + expectExportThrowsWithError(wsSameX, "two components have the observable 'nom_a'"); + } + { + RooGaussian tagged{"tagged", "tagged", nomA, a, RooFit::RooConst(1.0)}; + tagged.setAttribute("hs3_array_dist"); + RooWorkspace wsTagged; + wsTagged.import(tagged, RooFit::Silence()); + expectExportThrowsWithError(wsTagged, "is not a RooProdPdf"); + } +} + +// An array distribution paired with an unbinned data set with one entry (the form of the HS3 appendix on auxiliary +// measurements) gives a likelihood with one term per component. +TEST(RooFitHS3, ArrayDistLikelihood) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + const std::string json = R"({ + "metadata": {"hs3_version": "0.3"}, + "distributions": [ + {"name": "gauss_constraints", "type": "array_dist", "component_type": "gaussian_dist", + "default_parameters": {"sigma": 1.0}, + "components": [{"x": "nom_alpha_1", "mean": "alpha_1"}, {"x": "nom_alpha_2", "mean": "alpha_2"}, + {"x": "nom_alpha_3", "mean": "alpha_3", "sigma": 2.0}]} + ], + "data": [{"name": "global_observables", "type": "unbinned", "entries": [[0.5, -0.5, 1.0]], + "axes": [{"name": "nom_alpha_1", "min": -5.0, "max": 5.0}, {"name": "nom_alpha_2", "min": -5.0, "max": 5.0}, + {"name": "nom_alpha_3", "min": -5.0, "max": 5.0}]}], + "domains": [{"name": "default_domain", "type": "product_domain", "axes": [ + {"name": "nom_alpha_1", "min": -5.0, "max": 5.0}, {"name": "nom_alpha_2", "min": -5.0, "max": 5.0}, + {"name": "nom_alpha_3", "min": -5.0, "max": 5.0}, {"name": "alpha_1", "min": -5.0, "max": 5.0}, + {"name": "alpha_2", "min": -5.0, "max": 5.0}, {"name": "alpha_3", "min": -5.0, "max": 5.0}]}], + "parameter_points": [{"name": "default_values", "parameters": [ + {"name": "alpha_1", "value": 0.2}, {"name": "alpha_2", "value": 0.1}, {"name": "alpha_3", "value": -0.3}]}], + "likelihoods": [{"name": "likelihood", "distributions": ["gauss_constraints"], "data": ["global_observables"]}], + "analyses": [{"name": "analysis", "likelihood": "likelihood", "domains": ["default_domain"]}] + })"; + + RooWorkspace ws; + ASSERT_TRUE(RooJSONFactoryWSTool{ws}.importJSONfromString(json)); + auto *mc = dynamic_cast(ws.obj("analysis")); + ASSERT_NE(mc, nullptr); + ASSERT_NE(mc->GetPdf(), nullptr); + EXPECT_STREQ(mc->GetPdf()->GetName(), "gauss_constraints"); + auto *data = ws.data("global_observables"); + ASSERT_NE(data, nullptr); + ASSERT_EQ(data->numEntries(), 1); + + std::unique_ptr nll{mc->GetPdf()->createNLL(*data)}; + double expected = 0.0; + for (const char *name : {"nom_alpha_1_constraint", "nom_alpha_2_constraint", "nom_alpha_3_constraint"}) { + auto *component = ws.pdf(name); + ASSERT_NE(component, nullptr) << name; + RooArgSet observable; + component->getObservables(data->get(), observable); + expected -= std::log(component->getVal(observable)); + } + EXPECT_NEAR(nll->getVal(), expected, 1e-10 * std::abs(expected)); + + const std::string exported = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + RooWorkspace ws2; + ASSERT_TRUE(RooJSONFactoryWSTool{ws2}.importJSONfromString(exported)); + EXPECT_EQ(RooJSONFactoryWSTool{ws2}.exportJSONtoString(), exported); +} + +// A HistFactory modifier whose constraint is an array distribution gets the component that depends on its +// parameter. The model built from the array form equals the model built from one constraint per parameter, also +// for the constraint shared by both channels. +TEST(RooFitHS3, HistFactoryArrayDistConstraints) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + RooWorkspace individual; + ASSERT_TRUE(RooJSONFactoryWSTool{individual}.importJSONfromString(twoChannelHistFactoryJSON(false))); + RooWorkspace arrays; + const std::string warnings = captureMessages(RooFit::WARNING, RooFit::ObjectHandling, [&] { + ASSERT_TRUE(RooJSONFactoryWSTool{arrays}.importJSONfromString(twoChannelHistFactoryJSON(true))); + }); + EXPECT_TRUE(warnings.empty()) << warnings; + + // The channels hold the components, not the arrays. + auto constraintNames = [](RooWorkspace &ws, const char *channel) { + std::set names; + auto *model = dynamic_cast(ws.pdf(channel)); + if (model) { + for (RooAbsArg *pdf : model->pdfList()) { + if (!dynamic_cast(pdf)) { + names.insert(pdf->GetName()); + } + } + } + return names; + }; + EXPECT_EQ(constraintNames(arrays, "model_c1"), + (std::set{"alpha_a_constraint", "nom_gamma_c1_0_constraint", "nom_gamma_c1_1_constraint", + "nom_gs_0_constraint", "nom_gs_1_constraint"})); + EXPECT_EQ(constraintNames(arrays, "model_c2"), + (std::set{"alpha_a_constraint", "nom_gamma_c2_0_constraint", "nom_gamma_c2_1_constraint"})); + EXPECT_EQ(constraintNames(individual, "model_c1").size(), 5u); + EXPECT_EQ(constraintNames(individual, "model_c2").size(), 3u); + EXPECT_TRUE(arrays.pdf("pois_constraints")->getAttribute("hs3_array_dist")); + // The shared constraint is one object in both channels, so it is counted once. + auto countConstraints = [](RooWorkspace &ws) { + auto *mc = dynamic_cast(ws.obj("an")); + RooArgSet params; + mc->GetPdf()->getParameters(mc->GetObservables(), params); + return mc->GetPdf()->getAllConstraints(*mc->GetObservables(), params, true)->size(); + }; + EXPECT_EQ(countConstraints(arrays), 7u); + EXPECT_EQ(countConstraints(individual), 7u); + + const std::map point{ + {"mu", 1.3}, {"alpha_a", 0.7}, {"gamma_c1_0", 0.9}, {"gamma_c2_1", 1.1}, {"gs_0", 1.2}}; + const double nllIndividual = twoChannelNLL(individual, point); + const double nllArrays = twoChannelNLL(arrays, point); + EXPECT_FALSE(std::isnan(nllIndividual)); + EXPECT_NEAR(nllArrays, nllIndividual, 1e-10 * std::abs(nllIndividual)); +} + +// Errors and warnings when an array distribution does not fit the parameters that name it. +TEST(RooFitHS3, HistFactoryArrayDistConstraintErrors) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + auto document = [](std::string const &components) { + return R"({ + "metadata": {"hs3_version": "0.3"}, + "distributions": [ + {"name": "model_channel0", "type": "histfactory_dist", + "axes": [{"name": "x", "min": 0.0, "max": 2.0, "nbins": 2}], + "default_interpolation": {"type": "mult", "in": "poly6", "out": "exp"}, + "samples": [{"name": "sig", "data": {"contents": [10.0, 20.0]}, + "modifiers": [{"name": "a", "type": "normsys", "parameters": ["alpha_a"], + "data": {"lo": 0.9, "hi": 1.1}, "constraints": ["arr"]}]}]}, + {"name": "arr", "type": "array_dist", "component_type": "gaussian_dist", + "default_parameters": {"sigma": 1.0}, "components": [)" + + components + R"(]} + ], + "domains": [{"name": "default_domain", "type": "product_domain", "axes": [ + {"name": "alpha_a", "min": -5.0, "max": 5.0}, {"name": "b", "min": -5.0, "max": 5.0}, + {"name": "nom_1", "min": -5.0, "max": 5.0}, {"name": "nom_2", "min": -5.0, "max": 5.0}]}] + })"; + }; + + expectImportThrowsWithError(document(R"({"x": "nom_1", "mean": "b"})"), + "the array distribution 'arr' named as the constraint of 'alpha_a' (the normsys " + "modifier 'a' of sample 'sig') has no component that depends on this parameter"); + expectImportThrowsWithError(document(R"({"x": "nom_1", "mean": "alpha_a"}, {"x": "nom_2", "mean": "alpha_a"})"), + "has 2 components that depend on this parameter"); + + // A component that no parameter uses is reported: by the standard, the whole array enters the likelihood. + RooWorkspace ws; + const std::string warnings = captureMessages(RooFit::WARNING, RooFit::ObjectHandling, [&] { + EXPECT_TRUE(RooJSONFactoryWSTool{ws}.importJSONfromString( + document(R"({"x": "nom_1", "mean": "alpha_a"}, {"x": "nom_2", "mean": "b"})"))); + }); + EXPECT_NE(warnings.find("the components 'nom_2_constraint' of the array distribution 'arr' are not used"), + std::string::npos) + << warnings; + auto *model = dynamic_cast(ws.pdf("model_channel0")); + ASSERT_NE(model, nullptr); + EXPECT_NE(model->pdfList().find("nom_1_constraint"), nullptr); + EXPECT_EQ(model->pdfList().find("nom_2_constraint"), nullptr); +} + +// The HistFactory export collects the constraints of the modifiers into one array distribution per type and +// top-level pdf, which the modifiers name for each constrained parameter. The shared constraint is one component, the +// likelihood does not change under a round trip, and the second export equals the first. +TEST(RooFitHS3, HistFactoryArrayDistExport) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + RooWorkspace ws; + ASSERT_TRUE(RooJSONFactoryWSTool{ws}.importJSONfromString(twoChannelHistFactoryJSON(false))); + const std::string exported = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + + auto tree = RooFit::Detail::JSONTree::create(exported); + auto const &distributions = tree->rootnode()["distributions"]; + for (auto const &dist : distributions.children()) { + EXPECT_NE(dist["type"].val(), "gaussian_dist") << exported; + EXPECT_NE(dist["type"].val(), "poisson_dist") << exported; + } + const auto *gauss = RooJSONFactoryWSTool::findNamedChild(distributions, "an_simPdf_gaussian_constraints"); + ASSERT_NE(gauss, nullptr) << exported; + EXPECT_EQ((*gauss)["type"].val(), "array_dist"); + EXPECT_EQ((*gauss)["component_type"].val(), "gaussian_dist"); + ASSERT_EQ((*gauss)["components"].num_children(), 1u) << exported; + EXPECT_EQ((*gauss)["components"].child(0)["x"].val(), "alpha_a"); + const auto *pois = RooJSONFactoryWSTool::findNamedChild(distributions, "an_simPdf_poisson_constraints"); + ASSERT_NE(pois, nullptr) << exported; + EXPECT_EQ((*pois)["component_type"].val(), "poisson_dist"); + EXPECT_FALSE((*pois)["default_parameters"]["integer"].val_bool()); + // In the order of the channels, their samples and modifiers, and of the bins. + const std::vector gammas{"gs_0", "gs_1", "gamma_c1_0", "gamma_c1_1", "gamma_c2_0", "gamma_c2_1"}; + ASSERT_EQ((*pois)["components"].num_children(), gammas.size()) << exported; + for (std::size_t i = 0; i < gammas.size(); ++i) { + EXPECT_EQ((*pois)["components"].child(i)["x"].val(), "nom_" + gammas[i]) << exported; + EXPECT_EQ((*pois)["components"].child(i)["mean"].val(), gammas[i] + "_poisMean") << exported; + EXPECT_FALSE((*pois)["components"].child(i).has_child("integer")) << exported; + } + EXPECT_EQ(countOccurrences(exported, R"("constraints":["an_simPdf_gaussian_constraints"])"), 2u) << exported; + EXPECT_EQ( + countOccurrences(exported, R"("constraints":["an_simPdf_poisson_constraints","an_simPdf_poisson_constraints"])"), + 4u) + << exported; + + RooWorkspace ws2; + ASSERT_TRUE(RooJSONFactoryWSTool{ws2}.importJSONfromString(exported)); + const std::string exported2 = RooJSONFactoryWSTool{ws2}.exportJSONtoString(); + EXPECT_EQ(exported2, exported); + RooWorkspace ws3; + ASSERT_TRUE(RooJSONFactoryWSTool{ws3}.importJSONfromString(exported2)); + EXPECT_EQ(RooJSONFactoryWSTool{ws3}.exportJSONtoString(), exported); + + // The likelihood is unchanged by the round trip (evaluated last, since this moves the parameters). + const std::map point{ + {"mu", 1.3}, {"alpha_a", 0.7}, {"gamma_c1_0", 0.9}, {"gamma_c2_1", 1.1}, {"gs_0", 1.2}}; + const double nll = twoChannelNLL(ws, point); + EXPECT_FALSE(std::isnan(nll)); + EXPECT_NEAR(twoChannelNLL(ws2, point), nll, 1e-10 * std::abs(nll)); + + // Opt-out: one constraint distribution per parameter, no array_dist at all. + RooJSONFactoryWSTool::config().exportArrayDist = false; + const std::string plain = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + const std::string plainFromArrays = RooJSONFactoryWSTool{ws2}.exportJSONtoString(); + RooJSONFactoryWSTool::config().exportArrayDist = true; + EXPECT_EQ(plain.find("array_dist"), std::string::npos) << plain; + EXPECT_NE(plain.find(R"("constraints":["alpha_a_constraint"])"), std::string::npos) << plain; + EXPECT_NE(plain.find(R"("constraints":["gamma_c1_0_constraint","gamma_c1_1_constraint"])"), std::string::npos) + << plain; + RooWorkspace ws4; + ASSERT_TRUE(RooJSONFactoryWSTool{ws4}.importJSONfromString(plain)); + EXPECT_NEAR(twoChannelNLL(ws4, point), nll, 1e-10 * std::abs(nll)); + // From the workspace with imported arrays, the components are written on their own as well, and the arrays as + // product_dist. + EXPECT_EQ(plainFromArrays.find("array_dist"), std::string::npos) << plainFromArrays; + EXPECT_NE(plainFromArrays.find(R"("constraints":["nom_gamma_c1_0_constraint","nom_gamma_c1_1_constraint"])"), + std::string::npos) + << plainFromArrays; + RooWorkspace ws5; + ASSERT_TRUE(RooJSONFactoryWSTool{ws5}.importJSONfromString(plainFromArrays)); + EXPECT_NEAR(twoChannelNLL(ws5, point), nll, 1e-10 * std::abs(nll)); +} + +// The tutorial JSON imports and exports stably. Its Poisson constraints have a number as observable, so they cannot +// be components of an array distribution and are written on their own; the Gaussian constraints are collected. +TEST(RooFitHS3, TutorialRF515RoundTrip) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + const std::string path = std::string(gROOT->GetTutorialDir().Data()) + "/roofit/roofit/rf515_hfJSON.json"; + RooWorkspace ws; + ASSERT_TRUE(RooJSONFactoryWSTool{ws}.importJSON(path)); + const std::string exported = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + + auto tree = RooFit::Detail::JSONTree::create(exported); + auto const &distributions = tree->rootnode()["distributions"]; + const auto *poisson = RooJSONFactoryWSTool::findNamedChild(distributions, "gamma_bin_0_constraint"); + ASSERT_NE(poisson, nullptr) << exported; + EXPECT_EQ((*poisson)["type"].val(), "poisson_dist"); + // The top-level pdf is the RooSimultaneous "main" of the analysis. + const auto *gauss = RooJSONFactoryWSTool::findNamedChild(distributions, "main_gaussian_constraints"); + ASSERT_NE(gauss, nullptr) << exported; + EXPECT_EQ((*gauss)["type"].val(), "array_dist"); + EXPECT_EQ((*gauss)["components"].num_children(), 3u); + EXPECT_NE(exported.find(R"("constraints":["gamma_bin_0_constraint","gamma_bin_1_constraint"])"), std::string::npos) + << exported; + + RooWorkspace ws2; + ASSERT_TRUE(RooJSONFactoryWSTool{ws2}.importJSONfromString(exported)); + EXPECT_EQ(RooJSONFactoryWSTool{ws2}.exportJSONtoString(), exported); + + // None of the constraints has a global observable to pair with data (the Poissons have a number as observable, + // the Gaussians have the parameter as x and a number as mean), so all of them stay in "aux_distributions". + RooJSONFactoryWSTool::config().exportConstraintData = true; + const std::string withData = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + RooJSONFactoryWSTool::config().exportConstraintData = false; + auto withDataTree = RooFit::Detail::JSONTree::create(withData); + ASSERT_EQ(withDataTree->rootnode()["likelihoods"].num_children(), 1u) << withData; + const auto *likelihood = &withDataTree->rootnode()["likelihoods"].child(0); + ASSERT_EQ((*likelihood)["aux_distributions"].num_children(), 3u) << withData; + EXPECT_EQ((*likelihood)["aux_distributions"].child(0).val(), "main_gaussian_constraints") << withData; + EXPECT_EQ((*likelihood)["aux_distributions"].child(1).val(), "gamma_bin_0_constraint") << withData; + EXPECT_EQ((*likelihood)["aux_distributions"].child(2).val(), "gamma_bin_1_constraint") << withData; + EXPECT_EQ((*likelihood)["distributions"].num_children(), 1u) << withData; + EXPECT_EQ(withDataTree->rootnode()["data"].num_children(), 1u) << withData; +} + +// The constraint terms of the HistFactory channels are listed in the likelihood: in "aux_distributions" by default, +// with the observed values of their observables as constant parameters, or (config().exportConstraintData) in +// "distributions", each paired with a data set that holds the observed values. On import they enter the likelihood +// once, through the channels, and the observed values are taken from the data sets. +TEST(RooFitHS3, HistFactoryConstraintTermsInLikelihood) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + RooWorkspace ws; + ASSERT_TRUE(RooJSONFactoryWSTool{ws}.importJSONfromString(twoChannelHistFactoryJSON(false))); + // An observed value that differs from the default, to see it transported. + ws.var("nom_gamma_c1_0")->setVal(120.0); + const std::map point{ + {"mu", 1.3}, {"alpha_a", 0.7}, {"gamma_c1_0", 0.9}, {"gamma_c2_1", 1.1}, {"gs_0", 1.2}}; + + auto likelihood = [](std::string const &json) { + auto tree = RooFit::Detail::JSONTree::create(json); + const auto *node = RooJSONFactoryWSTool::findNamedChild(tree->rootnode()["likelihoods"], "an_likelihood"); + std::map> out; + if (node) { + for (const char *key : {"distributions", "data", "aux_distributions"}) { + if (node->has_child(key)) { + for (auto const &child : (*node)[key].children()) { + out[key].push_back(child.val()); + } + } + } + } + return out; + }; + using Names = std::vector; + + // Default: "aux_distributions". + const std::string exported = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + auto lh = likelihood(exported); + EXPECT_EQ(lh["distributions"], Names{"an_simPdf"}) << exported; + EXPECT_EQ(lh["data"], Names{"0"}) << exported; + EXPECT_EQ(lh["aux_distributions"], (Names{"an_simPdf_poisson_constraints", "an_simPdf_gaussian_constraints"})) + << exported; + { + RooWorkspace ws2; + ASSERT_TRUE(RooJSONFactoryWSTool{ws2}.importJSONfromString(exported)); + auto *mc = dynamic_cast(ws2.obj("an")); + ASSERT_NE(mc, nullptr); + // Not added again as external constraints: they are in the channels already. + EXPECT_TRUE(!mc->GetExternalConstraints() || mc->GetExternalConstraints()->empty()); + EXPECT_DOUBLE_EQ(ws2.var("nom_gamma_c1_0")->getVal(), 120.0); + EXPECT_TRUE(ws2.var("nom_gamma_c1_0")->isConstant()); + EXPECT_EQ(RooJSONFactoryWSTool{ws2}.exportJSONtoString(), exported); + const double nll = twoChannelNLL(ws, point); + EXPECT_NEAR(twoChannelNLL(ws2, point), nll, 1e-10 * std::abs(nll)); + } + + // Data route: the terms are paired with data sets holding the observed values. + RooJSONFactoryWSTool::config().exportConstraintData = true; + const std::string withData = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + RooJSONFactoryWSTool::config().exportConstraintData = false; + lh = likelihood(withData); + // The Gaussian of alpha_a has the parameter as x and a number as mean, so it has no global observable to pair + // with data and stays in "aux_distributions". + EXPECT_EQ(lh["distributions"], (Names{"an_simPdf", "an_simPdf_poisson_constraints"})) << withData; + EXPECT_EQ(lh["data"], (Names{"0", "an_simPdf_poisson_constraints_data"})) << withData; + EXPECT_EQ(lh["aux_distributions"], Names{"an_simPdf_gaussian_constraints"}) << withData; + { + auto tree = RooFit::Detail::JSONTree::create(withData); + const auto *data = + RooJSONFactoryWSTool::findNamedChild(tree->rootnode()["data"], "an_simPdf_poisson_constraints_data"); + ASSERT_NE(data, nullptr) << withData; + EXPECT_EQ((*data)["type"].val(), "unbinned"); + const std::vector gammas{"gs_0", "gs_1", "gamma_c1_0", "gamma_c1_1", "gamma_c2_0", "gamma_c2_1"}; + ASSERT_EQ((*data)["axes"].num_children(), gammas.size()) << withData; + ASSERT_EQ((*data)["entries"].num_children(), 1u) << withData; + ASSERT_EQ((*data)["entries"].child(0).num_children(), gammas.size()) << withData; + for (std::size_t i = 0; i < gammas.size(); ++i) { + auto const &axis = (*data)["axes"].child(i); + EXPECT_EQ(axis["name"].val(), "nom_" + gammas[i]) << withData; + EXPECT_DOUBLE_EQ(axis["min"].val_double(), 0.0); + EXPECT_DOUBLE_EQ(axis["max"].val_double(), 1000.0); + EXPECT_DOUBLE_EQ((*data)["entries"].child(0).child(i).val_double(), gammas[i] == "gamma_c1_0" ? 120.0 : 100.0); + } + EXPECT_EQ(RooJSONFactoryWSTool::findNamedChild(tree->rootnode()["data"], "an_simPdf_gaussian_constraints_data"), + nullptr) + << withData; + + RooWorkspace ws3; + ASSERT_TRUE(RooJSONFactoryWSTool{ws3}.importJSONfromString(withData)); + // The terms are not channels, and their data sets are not imported: the values live in the variables. + auto *mc = dynamic_cast(ws3.obj("an")); + ASSERT_NE(mc, nullptr); + auto *simPdf = dynamic_cast(mc->GetPdf()); + ASSERT_NE(simPdf, nullptr); + EXPECT_EQ(simPdf->indexCat().size(), 2); + EXPECT_EQ(ws3.data("an_simPdf_poisson_constraints_data"), nullptr); + EXPECT_DOUBLE_EQ(ws3.var("nom_gamma_c1_0")->getVal(), 120.0); + EXPECT_TRUE(ws3.var("nom_gamma_c1_0")->isConstant()); + ASSERT_NE(mc->GetGlobalObservables(), nullptr); + EXPECT_NE(mc->GetGlobalObservables()->find("nom_gamma_c1_0"), nullptr); + EXPECT_TRUE(!mc->GetExternalConstraints() || mc->GetExternalConstraints()->empty()); + const double nll = twoChannelNLL(ws, point); + EXPECT_NEAR(twoChannelNLL(ws3, point), nll, 1e-10 * std::abs(nll)); + RooJSONFactoryWSTool::config().exportConstraintData = true; + const std::string withData2 = RooJSONFactoryWSTool{ws3}.exportJSONtoString(); + RooJSONFactoryWSTool::config().exportConstraintData = false; + // The parameters were moved for the NLL, so compare the model and data sections. + for (const char *key : {"distributions", "functions", "data", "likelihoods"}) { + auto section = [key](std::string const &json) { + auto t = RooFit::Detail::JSONTree::create(json); + std::stringstream ss; + t->rootnode()[key].writeJSON(ss); + return ss.str(); + }; + EXPECT_EQ(section(withData2), section(withData)) << key; + } + } + + // Without arrays, every constraint is paired with its own data set. + RooJSONFactoryWSTool::config().exportArrayDist = false; + RooJSONFactoryWSTool::config().exportConstraintData = true; + const std::string plainWithData = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + RooJSONFactoryWSTool::config().exportArrayDist = true; + RooJSONFactoryWSTool::config().exportConstraintData = false; + lh = likelihood(plainWithData); + EXPECT_EQ(lh["distributions"].size(), 7u) << plainWithData; + EXPECT_EQ(lh["data"].size(), 7u) << plainWithData; + EXPECT_EQ(lh["distributions"][1], "gs_0_constraint") << plainWithData; + EXPECT_EQ(lh["data"][1], "gs_0_constraint_data") << plainWithData; + EXPECT_EQ(lh["aux_distributions"], Names{"alpha_a_constraint"}) << plainWithData; + { + auto tree = RooFit::Detail::JSONTree::create(plainWithData); + EXPECT_EQ(tree->rootnode()["data"].num_children(), 8u) << plainWithData; + RooWorkspace ws4; + ASSERT_TRUE(RooJSONFactoryWSTool{ws4}.importJSONfromString(plainWithData)); + const double nll = twoChannelNLL(ws, point); + EXPECT_NEAR(twoChannelNLL(ws4, point), nll, 1e-10 * std::abs(nll)); + } +} + +// Binned axes with an irregular binning carry "min" and "max" in addition to the "edges". +TEST(RooFitHS3, IrregularBinnedAxesHaveRange) +{ + RooRealVar x{"x", "x", 0.0, 3.0}; + std::vector edges{0.0, 1.0, 3.0}; + x.setBinning(RooBinning{2, edges.data()}); + RooDataHist dh{"dh", "dh", RooArgSet{x}}; + RooWorkspace ws; + ws.import(dh, RooFit::Silence()); + + const std::string json = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + auto tree = RooFit::Detail::JSONTree::create(json); + const auto *dataNode = RooJSONFactoryWSTool::findNamedChild(tree->rootnode()["data"], "dh"); + ASSERT_NE(dataNode, nullptr) << json; + const auto &axis = (*dataNode)["axes"].child(0); + EXPECT_TRUE(axis.has_child("edges")) << json; + EXPECT_DOUBLE_EQ(axis["min"].val_double(), 0.0) << json; + EXPECT_DOUBLE_EQ(axis["max"].val_double(), 3.0) << json; + + RooWorkspace ws2; + ASSERT_TRUE(RooJSONFactoryWSTool{ws2}.importJSONfromString(json)); + EXPECT_EQ(RooJSONFactoryWSTool{ws2}.exportJSONtoString(), json); +} + // Errors of zero in binned data stay zero, also before the first positive error. TEST(RooFitHS3, BinnedDataZeroErrors) { @@ -3534,6 +4699,151 @@ TEST(RooFitHS3, BinnedDataZeroErrors) EXPECT_DOUBLE_EQ(data->weightSquared(2), 0.0); } +// A HistFactory sample without "errors" contributes no statistical uncertainty to the staterror in HS3 >= 0.3. +// Before, RooFit assumed Poisson uncertainties for it. +TEST(RooFitHS3, HistFactorySampleWithoutErrors) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + auto makeJSON = [](std::string const &version) { + const std::string staterror = R"({"type": "staterror", "parameters": ["g0", "g1"], "constraint": "Gauss"})"; + return R"({ + "metadata": {"hs3_version": ")" + + version + R"("}, + "distributions": [ + { + "name": "model_channel0", + "type": "histfactory_dist", + "axes": [{"name": "x", "min": 0.0, "max": 2.0, "nbins": 2}], + "samples": [ + {"name": "a", "data": {"contents": [10.0, 20.0], "errors": [4.0, 6.0]}, "modifiers": [)" + + staterror + R"(]}, + {"name": "b", "data": {"contents": [10.0, 20.0]}, "modifiers": [)" + + staterror + R"(]} + ] + } + ] + })"; + }; + auto relativeErrors = [](RooWorkspace &ws) { + std::vector result; + for (const char *name : {"g0_sigma", "g1_sigma"}) { + auto *sigma = dynamic_cast(ws.arg(name)); + result.push_back(sigma ? sigma->getVal() : -1.0); + } + return result; + }; + + RooWorkspace ws; + ASSERT_TRUE(RooJSONFactoryWSTool{ws}.importJSONfromString(makeJSON("0.3"))); + const auto errors = relativeErrors(ws); + EXPECT_DOUBLE_EQ(errors[0], 4.0 / 20.0); + EXPECT_DOUBLE_EQ(errors[1], 6.0 / 40.0); + + RooWorkspace wsLegacy; + ASSERT_TRUE(RooJSONFactoryWSTool{wsLegacy}.importJSONfromString(makeJSON("0.2"))); + const auto legacyErrors = relativeErrors(wsLegacy); + EXPECT_DOUBLE_EQ(legacyErrors[0], std::sqrt(16.0 + 10.0) / 20.0); + EXPECT_DOUBLE_EQ(legacyErrors[1], std::sqrt(36.0 + 20.0) / 40.0); +} + +// HS3 has one per-bin modifier type, "shapefactor", whose parameters are constrained by the distributions named in +// "constraints", if given, and free otherwise. "shapesys" is a deprecated alias. HistFactory's ShapeSys (constrained) +// and ShapeFactor (free) are both written as "shapefactor". +TEST(RooFitHS3, HistFactoryShapeFactor) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + auto makeJSON = [](std::string const &version, std::string const &modifiers, + std::string const &constraintDist = "") { + return R"({ + "metadata": {"hs3_version": ")" + + version + R"("}, + "distributions": [ + { + "name": "model_channel0", + "type": "histfactory_dist", + "axes": [{"name": "x", "min": 0.0, "max": 2.0, "nbins": 2}], + "samples": [{"name": "sig", "data": {"contents": [10.0, 20.0]}, "modifiers": [)" + + modifiers + R"(]}] + })" + + constraintDist + R"( + ], + "domains": [{"name": "default_domain", "type": "product_domain", "axes": [ + {"name": "g_free_0", "min": 0, "max": 5}, {"name": "g_free_1", "min": 0, "max": 5}, + {"name": "g_c_0", "min": 0, "max": 5}, {"name": "g_c_1", "min": 0, "max": 5}]}] + })"; + }; + auto countGammaConstraints = [](RooWorkspace &ws) { + std::size_t n = 0; + for (RooAbsArg *pdf : ws.allPdfs()) { + if (dynamic_cast(pdf) || dynamic_cast(pdf)) { + ++n; + } + } + return n; + }; + const std::string freeModifier = + R"({"type": "shapefactor", "name": "free", "parameters": ["g_free_0", "g_free_1"]})"; + const std::string constrainedModifier = R"({"type": "shapefactor", "name": "constr", + "parameters": ["g_c_0", "g_c_1"], "constraints": ["g_c_0_constraint", null]})"; + const std::string json = + makeJSON("0.3", freeModifier + "," + constrainedModifier, + R"(, {"name": "g_c_0_constraint", "type": "gaussian_dist", "x": "g_c_0", "mean": 1.0, "sigma": 0.1})"); + + RooWorkspace ws; + ASSERT_TRUE(RooJSONFactoryWSTool{ws}.importJSONfromString(json)); + // The functions are named as in HistFactory: ShapeFactor for free, ShapeSys for constrained parameters. + EXPECT_NE(ws.function("channel0_free_shapeFactor"), nullptr); + EXPECT_NE(ws.function("channel0_constr_ShapeSys"), nullptr); + EXPECT_EQ(countGammaConstraints(ws), 1u); + ASSERT_NE(ws.var("g_free_0"), nullptr); + EXPECT_FALSE(ws.var("g_free_0")->isConstant()); + + // Both are written as "shapefactor", and free parameters don't prevent the round trip. + const std::string exported = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + EXPECT_EQ(countOccurrences(exported, R"("type":"shapefactor")"), 2u) << exported; + EXPECT_EQ(exported.find("shapesys"), std::string::npos) << exported; + EXPECT_NE(exported.find(R"("constraints":["model_channel0_gaussian_constraints",null])"), std::string::npos) + << exported; + RooWorkspace ws2; + ASSERT_TRUE(RooJSONFactoryWSTool{ws2}.importJSONfromString(exported)); + EXPECT_EQ(RooJSONFactoryWSTool{ws2}.exportJSONtoString(), exported); + + // "shapesys" is a deprecated alias. + std::string withAlias = json; + withAlias.replace(withAlias.find(R"("type": "shapefactor")"), std::string(R"("type": "shapefactor")").size(), + R"("type": "shapesys")"); + RooWorkspace ws3; + const std::string aliasWarnings = captureMessages(RooFit::WARNING, RooFit::ObjectHandling, [&] { + EXPECT_TRUE(RooJSONFactoryWSTool{ws3}.importJSONfromString(withAlias)); + }); + EXPECT_NE(aliasWarnings.find(R"(the modifier type "shapesys" is deprecated)"), std::string::npos) << aliasWarnings; + EXPECT_EQ(RooJSONFactoryWSTool{ws3}.exportJSONtoString(), exported); + + // In HS3 >= 0.3, "data.vals" without "constraints" don't create constraints. + RooWorkspace ws4; + const std::string valsWarnings = captureMessages(RooFit::WARNING, RooFit::ObjectHandling, [&] { + EXPECT_TRUE(RooJSONFactoryWSTool{ws4}.importJSONfromString( + makeJSON("0.3", R"({"type": "shapefactor", "name": "free", "parameters": ["g_free_0", "g_free_1"], + "data": {"vals": [0.1, 0.2]}})"))); + }); + EXPECT_NE(valsWarnings.find(R"(the "data.vals" of modifier 'free' of sample 'sig' are ignored)"), std::string::npos) + << valsWarnings; + EXPECT_NE(ws4.function("channel0_free_shapeFactor"), nullptr); + EXPECT_EQ(countGammaConstraints(ws4), 0u); + + // Before HS3 0.3, a "shapesys" with "data.vals" is constrained as before, without deprecation warning. + RooWorkspace ws5; + const std::string legacyWarnings = captureMessages(RooFit::WARNING, RooFit::ObjectHandling, [&] { + EXPECT_TRUE(RooJSONFactoryWSTool{ws5}.importJSONfromString( + makeJSON("0.2", R"({"type": "shapesys", "name": "constr", "data": {"vals": [0.1, 0.2]}})"))); + }); + EXPECT_EQ(legacyWarnings.find("deprecated"), std::string::npos) << legacyWarnings; + EXPECT_NE(ws5.function("channel0_constr_ShapeSys"), nullptr); + EXPECT_EQ(countGammaConstraints(ws5), 2u); +} + // A RooRealVar that is constant in all parameter points must stay a RooRealVar, and a RooConstVar must stay a // RooConstVar. In the parameter points both are written the same way (a value and "const": true). The only // difference is that the RooRealVar is listed in a domain and the RooConstVar is not: a parameter without domain is @@ -3783,11 +5093,16 @@ TEST(RooFitHS3, FormulaNormSysConstraint) RooWorkspace ws1; importFormulaNormSysChannel(ws1); const std::string json = RooJSONFactoryWSTool{ws1}.exportJSONtoString(); - EXPECT_NE(json.find(R"("constraint":"alpha_sysConstraint")"), std::string::npos) << json; + EXPECT_NE(json.find(R"("constraints":["model_ch_gaussian_constraints"])"), std::string::npos) << json; RooWorkspace ws2; ASSERT_TRUE(RooJSONFactoryWSTool{ws2}.importJSONfromString(json)); - EXPECT_NE(ws2.pdf("alpha_sysConstraint"), nullptr); + // The constraint survives the round trip; its name depends on the grouping of the constraints. + bool constrained = false; + for (RooAbsArg *pdf : ws2.allPdfs()) { + constrained |= dynamic_cast(pdf) && pdf->dependsOn(*ws2.var("alpha_sys")); + } + EXPECT_TRUE(constrained); } // The +1 and -1 sigma variations of a normsys that is encoded in a formula survive the round trip. @@ -4029,3 +5344,316 @@ TEST(RooFitHS3, WarningsAreShownByDefault) captureMessages(RooFit::WARNING, RooFit::ObjectHandling, [] { RooJSONFactoryWSTool::warning("test message"); }); EXPECT_NE(warnings.find("test message"), std::string::npos) << warnings; } + +// Binned contents are flattened with the last axis varying fastest (HS3 "Flattening of Binned Contents"). RooFit +// writes and reads them in that order. +TEST(RooFitHS3, FlatteningOrderOfBinnedContents) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + const std::string json = R"({ + "metadata": {"hs3_version": "0.3"}, + "data": [{"name": "d", "type": "binned", + "axes": [{"name": "x", "min": 0.0, "max": 2.0, "nbins": 2}, {"name": "y", "min": 0.0, "max": 3.0, "nbins": 3}], + "contents": [0.0, 1.0, 2.0, 3.0, 4.0, 5.0]}] + })"; + RooWorkspace ws; + ASSERT_TRUE(RooJSONFactoryWSTool{ws}.importJSONfromString(json)); + auto *data = dynamic_cast(ws.data("d")); + ASSERT_NE(data, nullptr); + RooRealVar &x = *ws.var("x"); + RooRealVar &y = *ws.var("y"); + for (int ix = 0; ix < 2; ++ix) { + for (int iy = 0; iy < 3; ++iy) { + x.setBin(ix); + y.setBin(iy); + EXPECT_DOUBLE_EQ(data->weight(RooArgSet{x, y}, 0), 3.0 * ix + iy) << "bin (" << ix << ", " << iy << ")"; + } + } + + // The export writes the same order back. + const std::string exported = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + EXPECT_NE(exported.find(R"("contents":[0,1,2,3,4,5])"), std::string::npos) << exported; +} + +// Categorical axes in domains create RooCategory objects, with one state per category. Binned and unbinned data can +// have categorical axes. Categories are identified by name: the data may list them in a different order than the +// domain, and the bins of binned data follow the order of the data axis. +TEST(RooFitHS3, CategoricalAxesInDomainsAndData) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + auto makeJSON = [](std::string const &data, std::string const &channelAxis = R"(["SR", "CR"])") { + return R"({ + "metadata": {"hs3_version": "0.3"}, + "domains": [{"name": "default_domain", "type": "product_domain", "axes": [ + {"name": "channel", "categories": )" + + channelAxis + R"(}, {"name": "x", "min": 0.0, "max": 2.0}]}], + "parameter_points": [{"name": "point", "parameters": [{"name": "channel", "value": "CR"}]}], + "data": [)" + + data + R"(] + })"; + }; + const std::string binned = R"({"name": "binned", "type": "binned", + "axes": [{"name": "channel", "categories": ["CR", "SR"]}, {"name": "x", "min": 0.0, "max": 2.0, "nbins": 2}], + "contents": [1.0, 2.0, 3.0, 4.0]})"; + const std::string unbinned = R"({"name": "unbinned", "type": "unbinned", + "axes": [{"name": "channel", "categories": ["SR", "CR"]}, {"name": "x", "min": 0.0, "max": 2.0}], + "entries": [["SR", 0.5], ["CR", 1.5], ["CR", 0.25]]})"; + + RooWorkspace ws; + ASSERT_TRUE(RooJSONFactoryWSTool{ws}.importJSONfromString(makeJSON(binned + "," + unbinned))); + + // The domain defines the category, with the position in "categories" as index. + RooCategory *channel = ws.cat("channel"); + ASSERT_NE(channel, nullptr); + EXPECT_EQ(channel->size(), 2u); + EXPECT_EQ(channel->lookupIndex("SR"), 0); + EXPECT_EQ(channel->lookupIndex("CR"), 1); + + // A parameter point sets the category by its label. + const RooArgSet *point = ws.getSnapshot("point"); + ASSERT_NE(point, nullptr); + auto *channelInPoint = dynamic_cast(point->find("channel")); + ASSERT_NE(channelInPoint, nullptr); + EXPECT_EQ(channelInPoint->getCurrentLabel(), std::string("CR")); + + // Binned data: the bins follow the order of the data axis ("CR" first), and are found by label. + auto *binnedData = dynamic_cast(ws.data("binned")); + ASSERT_NE(binnedData, nullptr); + RooRealVar &x = *ws.var("x"); + x.setBins(2); + const std::map, double> expected{ + {{"CR", 0}, 1.0}, {{"CR", 1}, 2.0}, {{"SR", 0}, 3.0}, {{"SR", 1}, 4.0}}; + for (auto const &item : expected) { + channel->setLabel(item.first.first); + x.setBin(item.first.second); + EXPECT_DOUBLE_EQ(binnedData->weight(RooArgSet{*channel, x}, 0), item.second) + << item.first.first << ", bin " << item.first.second; + } + + // Unbinned data: the entries of a categorical axis are the labels. + RooAbsData *unbinnedData = ws.data("unbinned"); + ASSERT_NE(unbinnedData, nullptr); + ASSERT_EQ(unbinnedData->numEntries(), 3); + const std::vector labels{"SR", "CR", "CR"}; + for (int i = 0; i < 3; ++i) { + auto const &row = *unbinnedData->get(i); + EXPECT_EQ(static_cast(row["channel"]).getCurrentLabel(), labels[i]) << "entry " << i; + } + + // Invalid documents. + expectImportThrowsWithError( + makeJSON(R"({"name": "u", "type": "unbinned", "axes": [{"name": "channel", "categories": ["SR", "CR"]}], + "entries": [["VR"]]})"), + "is not one of its categories"); + expectImportThrowsWithError( + makeJSON(R"({"name": "u", "type": "unbinned", "axes": [{"name": "channel", "categories": ["SR", "VR"]}], + "entries": [["SR"]]})"), + "differ from the states of the category"); + expectImportThrowsWithError( + makeJSON(R"({"name": "b", "type": "binned", "axes": [{"name": "channel", "categories": ["SR", "VR"]}], + "contents": [1.0, 2.0]})"), + "differ from those of the category"); + expectImportThrowsWithError(makeJSON("", R"(["SR", "SR"])"), "is listed twice"); + expectImportThrowsWithError(makeJSON("", R"([])"), "non-empty list of strings"); + expectImportThrowsWithError( + makeJSON(R"({"name": "b", "type": "binned", "axes": [{"name": "channel", "categories": ["SR", "CR"], "nbins": 2}], + "contents": [1.0, 2.0]})"), + "is categorical and must not have \"nbins\""); +} + +// RooFit writes per-bin arrays in the order HS3 prescribes for the axes it writes (last axis fastest): the contents of +// binned data and of RooHistFunc ("step"), and the parameters of ParamHistFunc ("step"). Each is built with its +// variables in both orders, and with a known value (or parameter name) per bin. +TEST(RooFitHS3, ExportedBinOrderFollowsAxes) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + RooRealVar x{"x", "x", 0.0, 2.0}; + x.setBins(2); + RooRealVar y{"y", "y", 0.0, 3.0}; + y.setBins(3); + auto valueOf = [](int ix, int iy) { return 10.0 * ix + iy + 1.0; }; + auto paramName = [](int ix, int iy) { return "g_" + std::to_string(ix) + "_" + std::to_string(iy); }; + + // Call `check(entry, ix, iy)` for the entry of `array` that HS3 assigns to the bin (ix, iy), given the axes of + // `node`. + auto forEachBin = [&](RooFit::Detail::JSONNode const &node, RooFit::Detail::JSONNode const &array, auto &&check) { + auto const &axes = node["axes"]; + ASSERT_EQ(axes.num_children(), 2u); + const bool xFirst = axes.child(0)["name"].val() == "x"; + EXPECT_EQ(axes.child(xFirst ? 1 : 0)["name"].val(), "y"); + ASSERT_EQ(array.num_children(), 6u); + for (int ix = 0; ix < 2; ++ix) { + for (int iy = 0; iy < 3; ++iy) { + const std::size_t k = xFirst ? ix * 3 + iy : iy * 2 + ix; + check(array.child(k), ix, iy); + } + } + }; + // Evaluate `name` in `ws` at the center of each bin. + auto expectValues = [&](RooWorkspace &ws, const char *name) { + RooAbsReal *func = ws.function(name); + ASSERT_NE(func, nullptr) << name; + for (int ix = 0; ix < 2; ++ix) { + for (int iy = 0; iy < 3; ++iy) { + ws.var("x")->setBin(ix); + ws.var("y")->setBin(iy); + EXPECT_DOUBLE_EQ(func->getVal(), valueOf(ix, iy)) << name << " bin (" << ix << ", " << iy << ")"; + } + } + }; + + for (bool xFirst : {true, false}) { + SCOPED_TRACE(xFirst ? "variables (x, y)" : "variables (y, x)"); + RooArgList vars = xFirst ? RooArgList{x, y} : RooArgList{y, x}; + + // Binned data and RooHistFunc. + RooDataHist dh{"dh", "dh", vars}; + for (int ix = 0; ix < 2; ++ix) { + for (int iy = 0; iy < 3; ++iy) { + x.setBin(ix); + y.setBin(iy); + dh.set(RooArgSet{x, y}, valueOf(ix, iy)); + } + } + RooHistFunc histFunc{"hist_func", "hist_func", vars, dh}; + + // ParamHistFunc: its parameters are ordered like the bins of a TH1, with the first variable varying fastest. + std::vector> params; + RooArgList paramList; + auto &first = static_cast(vars[0]); + auto &second = static_cast(vars[1]); + for (int i1 = 0; i1 < second.getBins(); ++i1) { + for (int i0 = 0; i0 < first.getBins(); ++i0) { + const int ix = xFirst ? i0 : i1; + const int iy = xFirst ? i1 : i0; + params.push_back(std::make_unique(paramName(ix, iy).c_str(), "", valueOf(ix, iy), 0.0, 100.0)); + paramList.add(*params.back()); + } + } + ParamHistFunc paramHistFunc{"param_hist_func", "param_hist_func", vars, paramList}; + // Check the assumption about ParamHistFunc that this test relies on. + for (int ix = 0; ix < 2; ++ix) { + for (int iy = 0; iy < 3; ++iy) { + x.setBin(ix); + y.setBin(iy); + ASSERT_DOUBLE_EQ(paramHistFunc.getVal(), valueOf(ix, iy)) << "bin (" << ix << ", " << iy << ")"; + } + } + + RooWorkspace ws; + ws.import(dh, RooFit::Silence()); + ws.import(histFunc, RooFit::RecycleConflictNodes(), RooFit::Silence()); + ws.import(paramHistFunc, RooFit::RecycleConflictNodes(), RooFit::Silence()); + const std::string json = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + auto tree = RooFit::Detail::JSONTree::create(json); + auto &root = tree->rootnode(); + + auto const *dataNode = RooJSONFactoryWSTool::findNamedChild(root["data"], "dh"); + ASSERT_NE(dataNode, nullptr) << json; + forEachBin(*dataNode, (*dataNode)["contents"], [&](auto const &entry, int ix, int iy) { + EXPECT_DOUBLE_EQ(entry.val_double(), valueOf(ix, iy)) << "data bin (" << ix << ", " << iy << ")"; + }); + + auto const *histNode = RooJSONFactoryWSTool::findNamedChild(root["functions"], "hist_func"); + ASSERT_NE(histNode, nullptr) << json; + forEachBin((*histNode)["data"], (*histNode)["data"]["contents"], [&](auto const &entry, int ix, int iy) { + EXPECT_DOUBLE_EQ(entry.val_double(), valueOf(ix, iy)) << "RooHistFunc bin (" << ix << ", " << iy << ")"; + }); + + auto const *paramNode = RooJSONFactoryWSTool::findNamedChild(root["functions"], "param_hist_func"); + ASSERT_NE(paramNode, nullptr) << json; + forEachBin(*paramNode, (*paramNode)["parameters"], [&](auto const &entry, int ix, int iy) { + EXPECT_EQ(entry.val(), paramName(ix, iy)) << "ParamHistFunc bin (" << ix << ", " << iy << ")"; + }); + + // Read back, the functions have the same value in every bin. + RooWorkspace imported; + ASSERT_TRUE(RooJSONFactoryWSTool{imported}.importJSONfromString(json)); + expectValues(imported, "hist_func"); + expectValues(imported, "param_hist_func"); + } +} + +// The per-bin parameters of a multi-dimensional HistFactory modifier follow the order of HS3 (last axis fastest), like +// the bin contents. Each gamma gets the relative uncertainty of its own bin, and the export keeps the order. +TEST(RooFitHS3, HistFactoryMultiDimensionalBinOrder) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + auto gammaName = [](int ix, int iy) { return "g_" + std::to_string(ix) + "_" + std::to_string(iy); }; + auto content = [](int ix, int iy) { return 100.0 * (1 + 3 * ix + iy); }; + auto relError = [](int ix, int iy) { return 0.01 * (1 + 3 * ix + iy); }; + auto makeJSON = [&](std::string const &version) { + std::string params; + std::string contents; + std::string errors; + for (int ix = 0; ix < 2; ++ix) { + for (int iy = 0; iy < 3; ++iy) { + const std::string sep = params.empty() ? "" : ", "; + params += sep + "\"" + gammaName(ix, iy) + "\""; + contents += sep + std::to_string(content(ix, iy)); + errors += sep + std::to_string(content(ix, iy) * relError(ix, iy)); + } + } + return R"({ + "metadata": {"hs3_version": ")" + + version + R"("}, + "distributions": [{ + "name": "model_channel0", "type": "histfactory_dist", + "axes": [{"name": "x", "min": 0.0, "max": 2.0, "nbins": 2}, {"name": "y", "min": 0.0, "max": 3.0, "nbins": 3}], + "samples": [{"name": "sig", "data": {"contents": [)" + + contents + R"(], "errors": [)" + errors + R"(]}, + "modifiers": [{"type": "staterror", "parameters": [)" + + params + R"(], "constraint": "Gauss"}]}] + }] + })"; + }; + + RooWorkspace ws; + ASSERT_TRUE(RooJSONFactoryWSTool{ws}.importJSONfromString(makeJSON("0.3"))); + auto *phf = dynamic_cast(ws.function("mc_stat_channel0")); + ASSERT_NE(phf, nullptr); + for (int ix = 0; ix < 2; ++ix) { + for (int iy = 0; iy < 3; ++iy) { + // The constraint of each gamma has the relative uncertainty of its bin. + auto *sigma = dynamic_cast(ws.arg((gammaName(ix, iy) + "_sigma").c_str())); + ASSERT_NE(sigma, nullptr) << gammaName(ix, iy); + EXPECT_NEAR(sigma->getVal(), relError(ix, iy), 1e-12) << gammaName(ix, iy); + // In its bin, the ParamHistFunc takes the value of that gamma. + ws.var(gammaName(ix, iy))->setVal(1.0 + 0.001 * (1 + 3 * ix + iy)); + } + } + for (int ix = 0; ix < 2; ++ix) { + for (int iy = 0; iy < 3; ++iy) { + ws.var("x")->setBin(ix); + ws.var("y")->setBin(iy); + EXPECT_DOUBLE_EQ(phf->getVal(), ws.var(gammaName(ix, iy))->getVal()) << "bin (" << ix << ", " << iy << ")"; + } + } + for (int ix = 0; ix < 2; ++ix) { + for (int iy = 0; iy < 3; ++iy) { + ws.var(gammaName(ix, iy))->setVal(1.0); + } + } + + // The export writes the gammas in the same order. + const std::string exported = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + EXPECT_NE(exported.find(R"("parameters":["g_0_0","g_0_1","g_0_2","g_1_0","g_1_1","g_1_2"])"), std::string::npos) + << exported; + + // Documents before HS3 0.3 don't define the order; they are rejected. + expectImportThrowsWithError(makeJSON("0.2"), "Documents before HS3 0.3 don't define the order"); + expectImportThrowsWithError(R"({ + "metadata": {"hs3_version": "0.2"}, + "domains": [{"name": "default_domain", "type": "product_domain", "axes": [ + {"name": "x", "min": 0.0, "max": 2.0}, {"name": "y", "min": 0.0, "max": 3.0}, + {"name": "a", "min": 0.0, "max": 5.0}, {"name": "b", "min": 0.0, "max": 5.0}]}], + "parameter_points": [{"name": "default_values", "parameters": [{"name": "a", "value": 1.0}, {"name": "b", "value": 2.0}]}], + "functions": [{"name": "phf", "type": "step", "variables": ["x", "y"], "parameters": ["a", "b"], + "axes": [{"name": "x", "min": 0.0, "max": 2.0, "nbins": 2}, {"name": "y", "min": 0.0, "max": 3.0, "nbins": 1}]}] + })", + "Documents before HS3 0.3 don't define the order"); +} diff --git a/tutorials/roofit/roofit/rf515_hfJSON.json b/tutorials/roofit/roofit/rf515_hfJSON.json index 382b8f9c39c06..2a65966ff8944 100644 --- a/tutorials/roofit/roofit/rf515_hfJSON.json +++ b/tutorials/roofit/roofit/rf515_hfJSON.json @@ -55,6 +55,39 @@ ] } ], + "domains": [ + { + "name": "staterror_parameters", + "type": "product_domain", + "axes": [ + { + "name": "gamma_bin_0", + "min": 0.0, + "max": 5.0 + }, + { + "name": "gamma_bin_1", + "min": 0.0, + "max": 5.0 + } + ] + } + ], + "parameter_points": [ + { + "name": "default_values", + "parameters": [ + { + "name": "gamma_bin_0", + "value": 1.0 + }, + { + "name": "gamma_bin_1", + "value": 1.0 + } + ] + } + ], "distributions": [ { "name": "model_channel1", @@ -85,12 +118,20 @@ "hi": 1.05, "lo": 0.95 }, - "constraint": "Gauss", + "constraints": [ + "alpha_syst1_constraint" + ], "name": "syst1", + "parameters": [ + "alpha_syst1" + ], "type": "normsys" }, { "name": "mu", + "parameters": [ + "mu" + ], "type": "normfactor" } ], @@ -113,12 +154,25 @@ "hi": 1.05, "lo": 0.95 }, - "constraint": "Gauss", + "constraints": [ + "alpha_syst2_constraint" + ], "name": "syst2", + "parameters": [ + "alpha_syst2" + ], "type": "normsys" }, { + "constraints": [ + "gamma_bin_0_constraint", + "gamma_bin_1_constraint" + ], "name": "mcstat", + "parameters": [ + "gamma_bin_0", + "gamma_bin_1" + ], "type": "staterror" } ], @@ -141,12 +195,25 @@ "hi": 1.05, "lo": 0.95 }, - "constraint": "Gauss", + "constraints": [ + "alpha_syst3_constraint" + ], "name": "syst3", + "parameters": [ + "alpha_syst3" + ], "type": "normsys" }, { + "constraints": [ + "gamma_bin_0_constraint", + "gamma_bin_1_constraint" + ], "name": "mcstat", + "parameters": [ + "gamma_bin_0", + "gamma_bin_1" + ], "type": "staterror" } ], @@ -154,6 +221,59 @@ } ], "type": "histfactory_dist" + }, + { + "name": "alpha_syst1_constraint", + "type": "gaussian_dist", + "x": "alpha_syst1", + "mean": 0.0, + "sigma": 1.0 + }, + { + "name": "alpha_syst2_constraint", + "type": "gaussian_dist", + "x": "alpha_syst2", + "mean": 0.0, + "sigma": 1.0 + }, + { + "name": "alpha_syst3_constraint", + "type": "gaussian_dist", + "x": "alpha_syst3", + "mean": 0.0, + "sigma": 1.0 + }, + { + "name": "gamma_bin_0_constraint", + "type": "poisson_dist", + "x": 400, + "mean": "gamma_bin_0_expected", + "integer": false + }, + { + "name": "gamma_bin_1_constraint", + "type": "poisson_dist", + "x": 100, + "mean": "gamma_bin_1_expected", + "integer": false + } + ], + "functions": [ + { + "name": "gamma_bin_0_expected", + "type": "product", + "factors": [ + "gamma_bin_0", + 400 + ] + }, + { + "name": "gamma_bin_1_expected", + "type": "product", + "factors": [ + "gamma_bin_1", + 100 + ] } ], "data": [ @@ -175,6 +295,6 @@ } ], "metadata": { - "hs3_version": "0.1.90" + "hs3_version": "0.3" } }