diff --git a/roofit/histfactory/test/testHistFactory.cxx b/roofit/histfactory/test/testHistFactory.cxx index 29acd2c4ded64..b8c6718e0da85 100644 --- a/roofit/histfactory/test/testHistFactory.cxx +++ b/roofit/histfactory/test/testHistFactory.cxx @@ -32,6 +32,7 @@ #include "../../roofitcore/test/gtest_wrapper.h" +#include #include #include @@ -700,15 +701,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 +741,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 +755,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 @@ -919,3 +916,77 @@ TEST(HistFactory, ShapeFactorValueAndRange) } } } + +/// The reference workspace of ROOT 6.16 is the standard case of a RooSimultaneous with combined data sets: it is +/// written as a "combined_dist" and "combined" data sets (HS3 #119), read back into the same objects, and the +/// likelihood is the same. The document is stable from the third export on: the 6.16 workspace has "binWidth_*" +/// constants that the imported model doesn't use, and the HistFactory importer sets the default parameter errors +/// (step widths), which the 6.16 workspace lacks. Both are unrelated to the combined objects. +TEST(HistFactory, HS3RoundTripOfReferenceWorkspace) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + std::unique_ptr file(TFile::Open("./ref_6.16_example_UsingC_combined_meas_model.root")); + ASSERT_TRUE(file && file->IsOpen()); + RooWorkspace *ws = nullptr; + file->GetObject("combined", ws); + ASSERT_NE(ws, nullptr); + auto *mc = dynamic_cast(ws->obj("ModelConfig")); + ASSERT_NE(mc, nullptr); + + const std::string json = RooJSONFactoryWSTool{*ws}.exportJSONtoString(); + EXPECT_NE(json.find(R"("type":"combined_dist")"), std::string::npos); + EXPECT_NE(json.find(R"("type":"combined")"), std::string::npos); + EXPECT_EQ(json.find("combined_datasets"), std::string::npos); + EXPECT_EQ(json.find("combined_distributions"), std::string::npos); + + RooWorkspace wsFromJson("new"); + RooJSONFactoryWSTool{wsFromJson}.importJSONfromString(json); + + auto *newmc = dynamic_cast(wsFromJson.obj("ModelConfig")); + ASSERT_NE(newmc, nullptr); + auto *newpdf = dynamic_cast(newmc->GetPdf()); + ASSERT_NE(newpdf, nullptr); + EXPECT_STREQ(newpdf->indexCat().GetName(), "channelCat"); + EXPECT_NE(newmc->GetObservables()->find("channelCat"), nullptr); + EXPECT_EQ(wsFromJson.allData().size(), ws->allData().size()); + for (const char *dataName : {"obsData", "asimovData"}) { + RooAbsData *data = ws->data(dataName); + RooAbsData *newdata = wsFromJson.data(dataName); + ASSERT_NE(data, nullptr) << dataName; + ASSERT_NE(newdata, nullptr) << dataName; + EXPECT_EQ(newdata->numEntries(), data->numEntries()) << dataName; + EXPECT_DOUBLE_EQ(newdata->sumEntries(), data->sumEntries()) << dataName; + EXPECT_NE(newdata->get()->find("channelCat"), nullptr) << dataName; + + // 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. + using namespace RooFit; + std::unique_ptr nll{mc->GetPdf()->createNLL(*data, GlobalObservables(*mc->GetGlobalObservables()))}; + std::unique_ptr newnll{ + newpdf->createNLL(*newdata, GlobalObservables(*newmc->GetGlobalObservables()))}; + const double nominal = nll->getVal(); + const double newNominal = newnll->getVal(); + const double alpha = ws->var("alpha_syst2")->getVal(); + const double mu = ws->var("SigXsecOverSM")->getVal(); + for (RooWorkspace *w : {ws, &wsFromJson}) { + w->var("alpha_syst2")->setVal(alpha + 0.5); + w->var("SigXsecOverSM")->setVal(mu + 0.5); + } + const double difference = nll->getVal() - nominal; + EXPECT_NE(difference, 0.0) << dataName; + EXPECT_NEAR(newnll->getVal() - newNominal, difference, 1e-8 * std::max(1.0, std::abs(difference))) << dataName; + for (RooWorkspace *w : {ws, &wsFromJson}) { + w->var("alpha_syst2")->setVal(alpha); + w->var("SigXsecOverSM")->setVal(mu); + } + } + + const std::string json2 = RooJSONFactoryWSTool{wsFromJson}.exportJSONtoString(); + RooWorkspace wsFromJson2("new2"); + RooJSONFactoryWSTool{wsFromJson2}.importJSONfromString(json2); + const std::string json3 = RooJSONFactoryWSTool{wsFromJson2}.exportJSONtoString(); + RooWorkspace wsFromJson3("new3"); + RooJSONFactoryWSTool{wsFromJson3}.importJSONfromString(json3); + EXPECT_EQ(RooJSONFactoryWSTool{wsFromJson3}.exportJSONtoString(), json3); +} diff --git a/roofit/hs3/inc/RooFitHS3/RooJSONFactoryWSTool.h b/roofit/hs3/inc/RooFitHS3/RooJSONFactoryWSTool.h index c711a16b355cf..9a4ff858797cf 100644 --- a/roofit/hs3/inc/RooFitHS3/RooJSONFactoryWSTool.h +++ b/roofit/hs3/inc/RooFitHS3/RooJSONFactoryWSTool.h @@ -21,6 +21,7 @@ #include #include +#include #include #include #include @@ -36,6 +37,9 @@ namespace RooStats { class ModelConfig; } class RooRealVar; +class RooAbsPdf; +class RooCategory; +class RooSimultaneous; class RooJSONFactoryWSTool { public: @@ -45,6 +49,13 @@ class RooJSONFactoryWSTool { bool allowExportInvalidNames = true; bool allowSanitizeNames = true; bool importNoDomainParametersAsRooConstVars = true; + /// HS3 requires the components of a combined distribution to have distinct variates, e.g. "mass_SR" and + /// "mass_CR" instead of one "mass". If enabled, the observables of every combined data set are renamed on + /// export to "obs_" (the labels of the key joined by "_"), or "obs__" if a component + /// has several observables; the components of the matching RooSimultaneous are written with the renamed + /// observables (functions depending on them get the suffix "_"). This allows exporting workspaces whose + /// components share an observable. Off by default: without it, such workspaces can't be exported. + bool renameCombinedObservables = false; }; static Config &config(); @@ -54,7 +65,7 @@ class RooJSONFactoryWSTool { struct CombinedData { std::string name; - std::map components; + std::vector columns; }; RooJSONFactoryWSTool(RooWorkspace &ws); @@ -184,14 +195,13 @@ class RooJSONFactoryWSTool { } static void exportAxis(RooFit::Detail::JSONNode &obsNode, RooRealVar const &var); + static void exportCategoricalAxis(RooFit::Detail::JSONNode &obsNode, RooCategory const &cat); static void exportHisto(RooArgSet const &vars, std::size_t n, double const *contents, RooFit::Detail::JSONNode &output); static void exportArray(std::size_t n, double const *contents, RooFit::Detail::JSONNode &output); - void exportCategory(RooAbsCategory const &cat, RooFit::Detail::JSONNode &node); - void queueExport(RooAbsArg const &arg) { _serversToExport.push_back(&arg); } void queueExportTemporary(RooAbsArg *arg) { @@ -228,6 +238,7 @@ class RooJSONFactoryWSTool { } } + void exportSimultaneous(RooSimultaneous const &simPdf, std::set &exportedObjectNames); void exportData(RooAbsData const &data); RooJSONFactoryWSTool::CombinedData exportCombinedData(RooAbsData const &data); @@ -246,8 +257,7 @@ class RooJSONFactoryWSTool { const std::vector &single); void exportSingleModelConfig(RooFit::Detail::JSONNode &rootnode, RooStats::ModelConfig const &mc, - std::string const &analysisName, - std::map const *dataComponents); + std::string const &analysisName, std::string const *dataName); // member variables const RooFit::Detail::JSONNode *_rootnodeInput = nullptr; @@ -261,6 +271,15 @@ class RooJSONFactoryWSTool { std::unique_ptr _domains; std::vector _serversToExport; std::vector _serversToDelete; + // The columns of the exported combined distributions: they are not parameters, and don't appear in parameter points. + std::set _columnCategories; + // Renaming of the observables of combined data sets, see Config::renameCombinedObservables: for each data set and + // each selection label, the new name of each observable; for each RooSimultaneous and selection label, the + // component with renamed observables to export instead of the original one. The objects are owned here. + std::map>> _combinedObservableRenames; + std::map> _renamedComponents; + std::vector> _renamedObjects; + void prepareCombinedObservableRenaming(); // 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..70004b4f31d4d 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); @@ -71,6 +120,17 @@ void Domains::readVariable(RooRealVar const &var) } } +void Domains::readCategory(RooAbsCategory const &cat) +{ + // The states in their insertion order: this is the order of the bins of a RooDataHist over the category, and the + // position becomes the index when the document is read back. + std::vector labels; + for (unsigned int i = 0; i < cat.size(); ++i) { + labels.push_back(cat.getOrdinal(i).first); + } + _map[defaultDomainName].readCategory(cat.GetName(), labels); +} + void Domains::writeVariable(RooRealVar &var) const { auto default_domain = _map.find(defaultDomainName); @@ -181,11 +241,29 @@ void Domains::ProductDomain::writeBinning(RooFit::Detail::JSONNode &node, Produc node["nbins"] << elem.nBins; } } +void Domains::ProductDomain::readCategory(const char *name, std::vector const &labels) +{ + auto &elem = _map[name]; + if (elem.hasMin || elem.hasMax || elem.hasNBins || !elem.edges.empty()) { + RooJSONFactoryWSTool::error(std::string("the variable '") + name + + "' is categorical, but the domain already lists it as continuous"); + } + if (!elem.categories.empty() && elem.categories != labels) { + RooJSONFactoryWSTool::error(std::string("the category '") + name + + "' is written to the domain twice, with different states"); + } + elem.categories = labels; +} + 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 +293,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 +352,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 +384,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..53654642d1387 100644 --- a/roofit/hs3/src/Domains.h +++ b/roofit/hs3/src/Domains.h @@ -18,6 +18,7 @@ #include class RooAbsBinning; +class RooAbsCategory; class RooRealVar; class RooWorkspace; @@ -35,12 +36,14 @@ namespace Detail { class Domains { public: void readVariable(RooRealVar const &); + void readCategory(RooAbsCategory const &); void writeVariable(RooRealVar &) const; void readJSON(RooFit::Detail::JSONNode const &); 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; @@ -48,12 +51,14 @@ class Domains { public: void readVariable(const RooRealVar &); void readVariable(const char *name, RooAbsBinning const &binning); + void readCategory(const char *name, std::vector const &labels); void writeVariable(RooRealVar &) const; void readJSON(RooFit::Detail::JSONNode const &); 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 +72,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..c96e99db5283f 100644 --- a/roofit/hs3/src/JSONFactories_HistFactory.cxx +++ b/roofit/hs3/src/JSONFactories_HistFactory.cxx @@ -444,6 +444,23 @@ void exportGammaParameters(JSONNode &mod, std::vector const ¶m mod["parameters"].fill_seq(paramNames); } +/// 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) +{ + if (std::none_of(constraints.begin(), constraints.end(), [](auto *pdf) { return pdf != nullptr; })) { + return; + } + auto &names = mod["constraints"].set_seq(); + for (auto *constraint : constraints) { + if (constraint) { + names.append_child() << constraint->GetName(); + } else { + names.append_child().set_null(); + } + } +} + RooRealVar &createNominal(RooWorkspace &ws, std::string const &parname, double val, double min, double max) { RooRealVar &nom = getOrCreate(ws, "nom_" + parname, val, min, max); @@ -486,8 +503,8 @@ RooAbsPdf &createLegacyConstraint(RooJSONFactoryWSTool &tool, const JSONNode &mo { 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 +531,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,9 +573,72 @@ 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) { + // HS3 references the constraint distributions by name. + if (auto name = singleConstraintName(mod); name && !name->empty()) { + auto *constraint = findNamedConstraint(tool, *name, sample); + 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(); @@ -600,6 +683,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 +704,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 +756,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 +796,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 +814,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 +840,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 +856,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()) { @@ -826,9 +914,26 @@ bool importHistSample(RooJSONFactoryWSTool &tool, RooDataHist &dh, RooArgSet con "not a supported legacy constraint type"); } } + // 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(); + } + const bool constrained = createConstraints ? !vals.empty() && constraint != "Const" + : std::any_of(constraintPdfs.begin(), constraintPdfs.end(), + [](RooAbsPdf *pdf) { return pdf != nullptr; }); + // Follow the naming of HistFactory: "__ShapeSys" for constrained and + // "__shapeFactor" for free parameters. + const std::string funcName = channelName + "_" + sysname + (constrained ? "_ShapeSys" : "_shapeFactor"); shapeElems.add(createPHF(funcName, sysname, parnames, vals, tool, constraints, varlist, constraint, - defaultGammaMin, defaultShapeSysGammaMax, minShapeUncertainty, - /*createConstraints=*/!hasConstraintList)); + defaultGammaMin, defaultShapeSysGammaMax, minShapeUncertainty, createConstraints)); for (auto *constraintPdf : constraintPdfs) { if (constraintPdf) { constraints.add(*constraintPdf); @@ -901,6 +1006,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 +1028,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 +1043,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 +1071,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 +1094,41 @@ 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"); + } + for (auto const &constraintName : *statErrConstraintNames) { + if (!constraintName) { + continue; + } + auto *constraint = findNamedConstraint(*tool, *constraintName, 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 +1306,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,6 +1469,7 @@ struct Sample { std::vector otherElements; bool useBarlowBeestonLight = false; std::vector staterrorParameters; + std::vector staterrorConstraints; Sample(const std::string &n) : name{n} {} }; @@ -1348,7 +1495,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 +1524,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(); @@ -1508,44 +1667,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 +1837,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) @@ -1753,12 +1903,10 @@ std::optional defaultInterpolation(const Channel &channel) bool exportChannel(RooJSONFactoryWSTool *tool, const Channel &channel, 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}); }; auto addModifier = [](JSONNode &modifiers, std::string const &name, const char *type) -> JSONNode & { auto &mod = modifiers.append_child(); @@ -1769,10 +1917,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 +1934,19 @@ 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 +1954,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 +1969,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); } for (const auto &other : sample.otherElements) { @@ -1854,9 +1985,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); } if (!observablesWritten) { @@ -2016,6 +2145,11 @@ bool tryExportHistFactory(RooJSONFactoryWSTool *tool, const std::string &pdfname } } } + for (auto *constraint : sample.staterrorConstraints) { + if (constraint) { + tool->queueExport(*constraint); + } + } } // Export all the custom modifiers diff --git a/roofit/hs3/src/JSONFactories_RooFitCore.cxx b/roofit/hs3/src/JSONFactories_RooFitCore.cxx index ce7102885bc66..ebbff62067e55 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 @@ -55,6 +57,9 @@ #include +#include +#include +#include #include "JSONIOUtils.h" #include @@ -66,7 +71,9 @@ #include #include #include +#include #include +#include #include #include @@ -409,6 +416,60 @@ bool importAddition(RooJSONFactoryWSTool *tool, const JSONNode &p) return true; } +// A "combined_dist" is a RooSimultaneous. Its index category is the single column, or a RooSuperCategory over the +// columns in the order of the document; the columns are categories declared by categorical axes of the domains. A +// "category_conditional_dist" is the same over its categorical parameters, tagged so that it is written back as +// such (RooFit evaluates a RooSimultaneous whose index category is not an observable as the selected component). +template +bool importSimultaneous(RooJSONFactoryWSTool *tool, const JSONNode &p) +{ + std::string const name(RooJSONFactoryWSTool::name(p)); + const char *type = Conditional ? "category_conditional_dist" : "combined_dist"; + std::string const context = "the " + std::string(type) + " '" + name + "'"; + // Whether a combined distribution is extended follows from its components; HS3 forbids declaring it. + if (!Conditional && p.has_child("extended")) { + RooJSONFactoryWSTool::error(context + " must not have \"extended\": whether it is extended follows from its " + + "components"); + } + std::vector columns = + readSelectionColumns(*tool->workspace(), p, Conditional ? "categorical_parameter" : "columns", context); + + auto const *componentsNode = p.find("components"); + if (!componentsNode || !componentsNode->is_seq() || componentsNode->num_children() == 0) { + RooJSONFactoryWSTool::error("\"components\" of " + context + " must be a non-empty list"); + } + std::map pdfMap; + for (JSONNode const &component : componentsNode->children()) { + if (!component.is_map() || !component.has_child("dist")) { + RooJSONFactoryWSTool::error("every component of " + context + " must have a \"key\" and a \"dist\""); + } + std::vector const key = readSelectionKey(component, columns, context); + std::string const label = selectionLabel(key); + if (pdfMap.find(label) != pdfMap.end()) { + RooJSONFactoryWSTool::error("the key " + selectionKeyToString(key) + " is used by two components of " + + context); + } + pdfMap[label] = tool->request(component["dist"].val(), name); + } + + RooSimultaneous *simPdf = nullptr; + if (columns.size() == 1) { + simPdf = &tool->wsImport(RooSimultaneous{name.c_str(), name.c_str(), pdfMap, *columns.front()}); + } else { + RooArgSet inputs; + for (RooCategory *column : columns) { + inputs.add(*column); + } + std::string const indexName = name + "_index"; + RooSuperCategory indexCat{indexName.c_str(), indexName.c_str(), inputs}; + simPdf = &tool->wsImport(RooSimultaneous{name.c_str(), name.c_str(), pdfMap, indexCat}); + } + if (Conditional) { + simPdf->setAttribute(categoryConditionalTag); + } + return true; +} + bool importAddPdf(RooJSONFactoryWSTool *tool, const JSONNode &p) { std::string name(RooJSONFactoryWSTool::name(p)); @@ -789,10 +850,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 +936,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) { @@ -1329,7 +1539,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,6 +1673,8 @@ 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")); @@ -1466,6 +1682,8 @@ STATIC_EXECUTE([]() { registerImporter("product_dist", false); registerImporter("sum", false); registerImporter("mixture_dist", false); + registerImporter>("combined_dist", false); + registerImporter>("category_conditional_dist", false); registerImporter("mixture_resolution_model", false); registerImporter("binsampling_dist", false); registerImporter("binsampling", false); @@ -1500,6 +1718,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); 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..2b2c2cd3dceaf 100644 --- a/roofit/hs3/src/JSONIOUtils.cxx +++ b/roofit/hs3/src/JSONIOUtils.cxx @@ -1,9 +1,18 @@ #include "JSONIOUtils.h" #include +#include +#include +#include +#include +#include +#include #include +#include #include +#include +#include #include #include @@ -167,3 +176,220 @@ 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."); + } +} + +namespace { + +void sortByName(std::vector &columns) +{ + std::sort(columns.begin(), columns.end(), + [](RooCategory const *l, RooCategory const *r) { return std::strcmp(l->GetName(), r->GetName()) < 0; }); +} + +} // namespace + +std::vector selectionColumns(RooAbsCategory const &indexCat) +{ + std::vector columns; + auto addColumn = [&](RooAbsArg const &arg) { + auto *cat = dynamic_cast(&arg); + if (!cat) { + RooJSONFactoryWSTool::error("the index category '" + std::string(arg.GetName()) + "' is a " + arg.ClassName() + + ", but only a RooCategory or a RooSuperCategory of RooCategories " + + "can be written as the columns of a combined distribution"); + } + columns.push_back(cat); + }; + if (auto *super = dynamic_cast(&indexCat)) { + for (RooAbsArg const *input : super->inputCatList()) { + addColumn(*input); + } + } else { + addColumn(indexCat); + } + sortByName(columns); + return columns; +} + +std::vector selectionColumns(RooAbsCollection const &vars) +{ + std::vector columns; + for (RooAbsArg const *arg : vars) { + if (auto *cat = dynamic_cast(arg)) { + columns.push_back(cat); + } + } + sortByName(columns); + return columns; +} + +std::vector> selectionKeys(std::vector const &columns) +{ + std::vector> keys{{}}; + for (RooCategory const *column : columns) { + std::vector> next; + for (auto const &key : keys) { + for (unsigned int i = 0; i < column->size(); ++i) { + next.push_back(key); + next.back().push_back(column->getOrdinal(i).first); + } + } + keys = std::move(next); + } + return keys; +} + +std::string selectionLabel(std::vector const &key) +{ + if (key.size() == 1) { + return key.front(); + } + std::string label = "{"; + for (std::size_t i = 0; i < key.size(); ++i) { + label += (i == 0 ? "" : ";") + key[i]; + } + return label + "}"; +} + +std::string selectionKeyToString(std::vector const &key) +{ + std::string out = "["; + for (std::size_t i = 0; i < key.size(); ++i) { + out += (i == 0 ? "\"" : ", \"") + key[i] + "\""; + } + return out + "]"; +} + +std::string selectionLabel(RooAbsCategory const &indexCat, std::vector const &columns, + std::vector const &key) +{ + auto *super = dynamic_cast(&indexCat); + if (!super) { + return key.front(); + } + std::vector inInputOrder; + for (RooAbsArg const *input : super->inputCatList()) { + auto found = std::find_if(columns.begin(), columns.end(), + [&](RooCategory const *c) { return std::strcmp(c->GetName(), input->GetName()) == 0; }); + inInputOrder.push_back(key[found - columns.begin()]); + } + return selectionLabel(inInputOrder); +} + +std::vector +readSelectionColumns(RooWorkspace const &ws, JSONNode const &node, const char *field, std::string const &context) +{ + auto const *list = node.find(field); + if (!list || !list->is_seq() || list->num_children() == 0) { + RooJSONFactoryWSTool::error("\"" + std::string(field) + "\" of " + context + + " must be a non-empty list of names of categorical axes"); + } + std::vector columns; + for (JSONNode const &child : list->children()) { + if (child.is_container() || child.val().empty()) { + RooJSONFactoryWSTool::error("\"" + std::string(field) + "\" of " + context + + " must be a non-empty list of names of categorical axes"); + } + RooCategory *cat = ws.cat(child.val()); + if (!cat) { + RooJSONFactoryWSTool::error("the column '" + child.val() + "' of " + context + + " is not declared by a categorical axis of a domain"); + } + if (std::find(columns.begin(), columns.end(), cat) != columns.end()) { + RooJSONFactoryWSTool::error("the column '" + child.val() + "' of " + context + " is listed twice"); + } + columns.push_back(cat); + } + return columns; +} + +std::vector +readSelectionKey(JSONNode const &component, std::vector const &columns, std::string const &context) +{ + auto const *keyNode = component.find("key"); + if (!keyNode || !keyNode->is_seq() || keyNode->num_children() != columns.size()) { + RooJSONFactoryWSTool::error("every component of " + context + " must have a \"key\" with " + + std::to_string(columns.size()) + " labels, one per column"); + } + std::vector key; + std::size_t i = 0; + for (JSONNode const &label : keyNode->children()) { + if (label.is_container() || !columns[i]->hasLabel(label.val())) { + RooJSONFactoryWSTool::error("the label '" + (label.is_container() ? std::string("") : label.val()) + + "' in position " + std::to_string(i) + " of a key of " + context + + " is not one of the categories of '" + columns[i]->GetName() + "'"); + } + key.push_back(label.val()); + ++i; + } + return key; +} diff --git a/roofit/hs3/src/JSONIOUtils.h b/roofit/hs3/src/JSONIOUtils.h index 2ac14662be7dc..2365e766b605e 100644 --- a/roofit/hs3/src/JSONIOUtils.h +++ b/roofit/hs3/src/JSONIOUtils.h @@ -6,6 +6,13 @@ #include class RooAbsBinning; +class RooAbsCategory; +class RooAbsCollection; +class RooAbsReal; +class RooArgList; +class RooCategory; +class RooWorkspace; +class ParamHistFunc; bool startsWith(std::string_view str, std::string_view prefix); bool endsWith(std::string_view str, std::string_view suffix); @@ -17,4 +24,45 @@ 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); + +// A RooSimultaneous with this attribute is a "category_conditional_dist": its index category is a parameter, and it +// evaluates to the component selected by the current state ("switch mode"). Without it, a RooSimultaneous is a +// "combined_dist", whose index category is the column of the paired combined data set. RooFit decides between the +// two by the normalization set, i.e. by context; the attribute makes the export context-free and is set when a +// "category_conditional_dist" is imported. +inline constexpr const char *categoryConditionalTag = "hs3_category_conditional_dist"; + +// Selection keys (HS3 "combined" data sets and "combined_dist"): a key is one label per column, and the columns are +// categories declared by categorical axes of the domains. In RooFit, the columns are the index category of a +// RooSimultaneous (a RooCategory), or the input categories of a RooSuperCategory index, or the category observables +// of a combined data set. They are sorted by name, so that a distribution and a data set agree on the order of the +// columns without looking at each other. +std::vector selectionColumns(RooAbsCategory const &indexCat); +std::vector selectionColumns(RooAbsCollection const &vars); +// All keys of the columns, in the order of the columns and of their states, with the first column varying slowest. +std::vector> selectionKeys(std::vector const &columns); +// The state label of a RooCategory (one column) or of a RooSuperCategory over the columns in that order. +std::string selectionLabel(std::vector const &key); +// The state label of `indexCat` for a key with the labels in the order of `columns`: the label itself for a +// RooCategory, or "{a;b}" with the labels in the order of the inputs for a RooSuperCategory. +std::string selectionLabel(RooAbsCategory const &indexCat, std::vector const &columns, + std::vector const &key); +// A key as written in the document, for messages: ["a", "b"]. +std::string selectionKeyToString(std::vector const &key); +// Read the columns of a combined distribution or data set from the list `field` of `node`: each entry must be the +// name of a category in the workspace, i.e. declared by a categorical axis of a domain. +std::vector readSelectionColumns(RooWorkspace const &ws, RooFit::Detail::JSONNode const &node, + const char *field, std::string const &context); +// Read and validate the "key" of a component: one label per column, each one of the categories of its column. +std::vector readSelectionKey(RooFit::Detail::JSONNode const &component, + std::vector const &columns, 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..e0a45f9165fed 100644 --- a/roofit/hs3/src/RooJSONFactoryWSTool.cxx +++ b/roofit/hs3/src/RooJSONFactoryWSTool.cxx @@ -20,10 +20,13 @@ #include #include #include +#include #include #include #include +#include #include +#include #include #include #include @@ -34,6 +37,8 @@ #include #include +#include +#include #include #include #include @@ -106,7 +111,7 @@ tool.writedoc("hs3.tex") ~~~ */ -constexpr auto hs3VersionTag = "0.2"; +constexpr auto hs3VersionTag = "0.3"; using RooFit::Detail::JSONNode; using RooFit::Detail::JSONTree; @@ -123,12 +128,6 @@ std::vector valsToStringVec(JSONNode const &node) return out; } -// True if the number of components in `data` matches the number of categories in `pdf`. -bool matches(const RooJSONFactoryWSTool::CombinedData &data, const RooSimultaneous *pdf) -{ - return data.components.size() == pdf->indexCat().size(); -} - // True if the entire string parses as a number (integer or floating-point). bool isNumber(const std::string &str) { @@ -367,8 +366,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); @@ -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) { @@ -530,10 +595,78 @@ loadData(const JSONNode &p, RooWorkspace &workspace, RooFit::JSONIO::Detail::Dom RooJSONFactoryWSTool::error("RooJSONFactoryWSTool() failed to create dataset " + name); } +// What the document says about a combined data set: its columns and the keys with a sub-data set, each key as a +// map from the column to the label. Used to check the pairing with a combined distribution in a likelihood. +struct CombinedDataInfo { + std::vector columns; + std::set> keys; +}; + +// Check that a combined data set paired with a RooSimultaneous in a likelihood has a sub-data set for every component. +// RooFit has no notion of an unobserved component: a key without sub-data set is a channel without entries, which +// adds the expected number of events of an extended component to the likelihood. That changes the fit, so it is an +// error (https://github.com/root-project/root/issues/24004). For a non-extended component it changes nothing, and a +// sub-data set without component is ignored by RooFit; both are reported as warnings. +void checkCombinedPairing(RooSimultaneous const &simPdf, std::string const &dataName, CombinedDataInfo const &info) +{ + std::vector columns = selectionColumns(simPdf.indexCat()); + std::vector columnNames; + for (RooCategory const *column : columns) { + columnNames.push_back(column->GetName()); + } + std::vector dataColumns = info.columns; + std::sort(dataColumns.begin(), dataColumns.end()); + if (dataColumns != columnNames) { + RooJSONFactoryWSTool::error("the likelihood pairs '" + std::string(simPdf.GetName()) + + "' with the combined data set '" + dataName + "', but their columns differ"); + } + auto keyToString = [&](std::map const &key) { + std::vector labels; + for (std::string const &column : info.columns) { + labels.push_back(key.at(column)); + } + return selectionKeyToString(labels); + }; + std::set> pdfKeys; + for (std::vector const &labels : selectionKeys(columns)) { + RooAbsPdf const *component = simPdf.getPdf(selectionLabel(simPdf.indexCat(), columns, labels)); + if (!component) { + continue; + } + std::map key; + for (std::size_t i = 0; i < columns.size(); ++i) { + key[columnNames[i]] = labels[i]; + } + pdfKeys.insert(key); + if (info.keys.count(key) > 0) { + continue; + } + if (component->canBeExtended()) { + RooJSONFactoryWSTool::error( + "the combined data set '" + dataName + "' has no sub-data set for the key " + keyToString(key) + + " of the extended component '" + component->GetName() + "' of '" + simPdf.GetName() + + "'. RooFit cannot leave a component unobserved: it would add the expected number of events of the " + + "component to the likelihood, see https://github.com/root-project/root/issues/24004"); + } + RooJSONFactoryWSTool::warning("RooFitHS3: the combined data set '" + dataName + + "' has no sub-data set for the key " + keyToString(key) + " of the component '" + + component->GetName() + "' of '" + simPdf.GetName() + + "', which RooFit treats as present and empty."); + } + for (auto const &key : info.keys) { + if (pdfKeys.count(key) == 0) { + RooJSONFactoryWSTool::warning("RooFitHS3: the sub-data set for the key " + keyToString(key) + + " of the combined data set '" + dataName + "' is ignored, because '" + + simPdf.GetName() + "' has no component for it."); + } + } +} + // 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) + const std::vector> &datasets, + std::map const &combinedData) { // if this is a toplevel pdf, also create a modelConfig for it std::string const &analysisName = RooJSONFactoryWSTool::name(analysisNode); @@ -569,7 +702,8 @@ void importAnalysis(const JSONNode &rootnode, const JSONNode &analysisNode, cons for (auto &nameNode : (*nllNode)["data"].children()) { bool found = false; for (const auto &d : datasets) { - if (d->GetName() == nameNode.val()) { + // Components of combined data sets have been moved out of the list. + if (d && d->GetName() == nameNode.val()) { found = true; observables.add(*d->get(), true); } @@ -578,6 +712,16 @@ void importAnalysis(const JSONNode &rootnode, const JSONNode &analysisNode, cons throw std::runtime_error("dataset '" + nameNode.val() + "' cannot be found!"); } + // A combined distribution paired with a combined data set must have a sub-data set for every component. + std::vector const nllDataNames = valsToStringVec((*nllNode)["data"]); + for (std::size_t i = 0; i < std::min(nllDistNames.size(), nllDataNames.size()); ++i) { + auto *simPdf = dynamic_cast(workspace.pdf(nllDistNames[i])); + auto info = combinedData.find(nllDataNames[i]); + if (simPdf && !simPdf->getAttribute(categoryConditionalTag) && info != combinedData.end()) { + checkCombinedPairing(*simPdf, nllDataNames[i], info->second); + } + } + JSONNode const *pdfNameNode = mcAuxNode ? mcAuxNode->find("pdfName") : nullptr; std::string const pdfName = pdfNameNode ? pdfNameNode->val() : "simPdf"; @@ -609,15 +753,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; } - return out; + 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 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 +793,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()) { @@ -772,6 +936,151 @@ void sortByName(T &coll) std::sort(coll.begin(), coll.end(), [](auto &l, auto &r) { return strcmp(l->GetName(), r->GetName()) < 0; }); } +// Build the combined data sets of the document ("type": "combined"): a RooDataSet with one category column per entry +// of "columns" and the variables of all components, filled row by row from the components, which are then no longer +// imported on their own. The columns are categories declared by categorical axes of the domains. A key without a +// sub-data set is a state without entries: RooFit has no notion of an absent component. +// The names of the variables that appear in more than one of `datasets`. +template +std::vector sharedVariableNames(DataSets const &datasets, GetData getData) +{ + std::map counts; + for (auto const &item : datasets) { + for (RooAbsArg const *var : *getData(item)->get()) { + ++counts[var->GetName()]; + } + } + std::vector shared; + for (auto const &[name, count] : counts) { + if (count > 1) { + shared.push_back(name); + } + } + return shared; +} + +std::string joinQuoted(std::vector const &names) +{ + std::string out; + for (std::string const &name : names) { + out += (out.empty() ? "'" : ", '") + name + "'"; + } + return out; +} + +std::map assembleCombinedData(JSONNode const &dataNode, RooWorkspace const &workspace, + std::vector> &datasets) +{ + std::map infos; + for (JSONNode const &p : dataNode.children()) { + if (!p.is_map() || !p.has_child("type") || p["type"].val() != "combined") { + continue; + } + std::string const name = RooJSONFactoryWSTool::name(p); + RooJSONFactoryWSTool::testValidName(name, true); + std::string const context = "the combined data set '" + name + "'"; + if (p.has_child("axes")) { + RooJSONFactoryWSTool::error(context + " must not have \"axes\": its columns are declared by the domains"); + } + std::vector columns = readSelectionColumns(workspace, p, "columns", context); + auto const *componentsNode = p.find("components"); + if (!componentsNode || !componentsNode->is_seq() || componentsNode->num_children() == 0) { + RooJSONFactoryWSTool::error("\"components\" of " + context + " must be a non-empty list"); + } + + struct Component { + std::vector key; + std::unique_ptr data; + }; + std::vector components; + std::set labels; + for (JSONNode const &c : componentsNode->children()) { + if (!c.is_map() || !c.has_child("dataset")) { + RooJSONFactoryWSTool::error("every component of " + context + " must have a \"key\" and a \"dataset\""); + } + std::vector key = readSelectionKey(c, columns, context); + if (!labels.insert(selectionLabel(key)).second) { + RooJSONFactoryWSTool::error("the key " + selectionKeyToString(key) + " is used by two components of " + + context); + } + std::string const dataName = c["dataset"].val(); + auto found = + std::find_if(datasets.begin(), datasets.end(), [&](auto &d) { return d && d->GetName() == dataName; }); + if (found == datasets.end()) { + RooJSONFactoryWSTool::error("the component '" + dataName + "' of " + context + + " is not a data set of the document, or it is already a component of " + + "another combined data set"); + } + components.push_back({std::move(key), std::move(*found)}); + } + // HS3 requires distinct axis names, since the components of a combined distribution share no variate. RooFit + // can still read such a document. + std::vector shared = sharedVariableNames(components, [](auto const &c) { return c.data.get(); }); + if (!shared.empty()) { + RooJSONFactoryWSTool::warning("RooFitHS3: the sub-data sets of " + context + " share the axes " + + joinQuoted(shared) + + ". HS3 requires distinct axis names, because the components of a combined " + "distribution share no variate. RooFit reads the document anyway."); + } + + // The variables: the columns first, then those of the components. A weight variable if any component is + // weighted, with errors if any component is binned, like RooDataSet does when importing a map of data sets. + RooArgSet columnClones; // owns the clones; `vars` must not own anything, it also lists the components' variables + for (RooCategory const *column : columns) { + columnClones.addClone(*column); + } + RooArgSet vars{columnClones}; + bool weighted = false; + bool storeErrors = false; + std::string weightName = RooFit::WeightVar().getString(0); + for (auto const &c : components) { + vars.add(*c.data->get(), /*silent=*/true); + weighted |= c.data->isWeighted(); + if (auto *ds = dynamic_cast(c.data.get()); ds && ds->isWeighted()) { + weightName = ds->weightVar()->GetName(); + } + storeErrors |= dynamic_cast(c.data.get()) != nullptr; + } + std::unique_ptr combined; + if (weighted) { + RooRealVar weightVar{weightName.c_str(), weightName.c_str(), 1.0}; + vars.add(weightVar); + RooArgSet errorSet{weightVar}; + combined = std::make_unique(name.c_str(), name.c_str(), vars, RooFit::WeightVar(weightVar), + storeErrors ? RooFit::StoreError(errorSet) : RooCmdArg::none()); + } else { + combined = std::make_unique(name.c_str(), name.c_str(), vars); + } + for (auto const &c : components) { + for (std::size_t i = 0; i < columns.size(); ++i) { + static_cast(vars[columns[i]->GetName()]).setLabel(c.key[i].c_str()); + } + for (int i = 0; i < c.data->numEntries(); ++i) { + vars.assign(*c.data->get(i)); + if (weighted) { + combined->add(vars, c.data->weight(), storeErrors ? c.data->weightError(RooAbsData::SumW2) : 0.0); + } else { + combined->add(vars); + } + } + } + datasets.emplace_back(std::move(combined)); + + CombinedDataInfo &info = infos[name]; + for (RooCategory const *column : columns) { + info.columns.push_back(column->GetName()); + } + for (auto const &c : components) { + std::map key; + for (std::size_t i = 0; i < columns.size(); ++i) { + key[columns[i]->GetName()] = c.key[i]; + } + info.keys.insert(key); + } + } + return infos; +} + /// Check whether an object is a top-level object of the workspace, i.e. not a /// component of some other workspace object. Only clients that are part of the /// workspace count: RooFit also registers evaluation artifacts as clients, @@ -789,22 +1098,6 @@ bool isTopLevel(RooAbsArg const &arg, RooWorkspace const &ws) return true; } -/// Find the single category observable of a dataset, if any. -RooAbsCategory *findCategoryObservable(RooAbsData const &data) -{ - RooAbsCategory *cat = nullptr; - for (RooAbsArg *obs : *data.get()) { - if (auto *c = dynamic_cast(obs)) { - if (cat) { - RooJSONFactoryWSTool::error("dataset '" + std::string(data.GetName()) + - " has several category observables!"); - } - cat = c; - } - } - return cat; -} - /// Warn if the HS3 version of an imported document is not understood or newer than what RooFit writes. void checkHS3Version(std::string const &version) { @@ -825,8 +1118,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 +1145,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; @@ -1022,6 +1321,16 @@ RooAbsReal *RooJSONFactoryWSTool::requestImpl(const std::string &obj // Export a single variable (RooRealVar or RooConstVar) `v` as a named child of `node`. void RooJSONFactoryWSTool::exportVariable(const RooAbsArg *v, JSONNode &node, bool storeConstant, bool storeBins) { + // A categorical parameter: its value is the label. The columns of combined distributions are not parameters. + if (auto *cat = dynamic_cast(v)) { + if (_columnCategories.count(cat->GetName()) > 0) { + return; + } + JSONNode &var = appendNamedChild(node, cat->GetName()); + var["value"] << cat->getCurrentLabel(); + _domains->readCategory(*cat); + return; + } auto *cv = dynamic_cast(v); auto *rrv = dynamic_cast(v); if (!cv && !rrv) @@ -1090,27 +1399,11 @@ void RooJSONFactoryWSTool::exportObject(RooAbsArg const &func, std::set(&func)) { - // RooSimultaneous is not used in the HS3 standard, we only export the - // dependents and some ROOT internal information. - exportObjects(func.servers(), exportedObjectNames); - - std::vector channelNames; - for (auto const &item : simPdf->indexCat()) { - channelNames.push_back(item.first); - } - - auto &infoNode = getRooFitInternal(*_rootnodeOutput, "combined_distributions").set_map(); - auto &child = infoNode[simPdf->GetName()].set_map(); - child["index_cat"] << simPdf->indexCat().GetName(); - exportCategory(simPdf->indexCat(), child); - child["distributions"].set_seq(); - for (auto const &item : simPdf->indexCat()) { - child["distributions"].append_child() << simPdf->getPdf(item.first.c_str())->GetName(); - } - + exportSimultaneous(*simPdf, exportedObjectNames); return; } else if (dynamic_cast(&func)) { - // categories are created by the respective RooSimultaneous, so we're skipping the export here + // Categories are not objects of the document: the index category of a RooSimultaneous is declared by the + // categorical axes of the domain that its "columns" name. return; } else if (dynamic_cast(&func) || dynamic_cast(&func)) { exportVariable(&func, *_varsNode, true, false); @@ -1220,6 +1513,70 @@ void RooJSONFactoryWSTool::exportObject(RooAbsArg const &func, std::set &exportedObjectNames) +{ + RooAbsCategory const &indexCat = simPdf.indexCat(); + std::vector columns = selectionColumns(indexCat); + const bool conditional = simPdf.getAttribute(categoryConditionalTag); + + JSONNode &elem = appendNamedChild((*_rootnodeOutput)["distributions"], simPdf.GetName()); + elem["type"] << (conditional ? "category_conditional_dist" : "combined_dist"); + auto &columnsNode = elem[conditional ? "categorical_parameter" : "columns"].set_seq(); + for (RooCategory const *column : columns) { + columnsNode.append_child() << column->GetName(); + _domains->readCategory(*column); + if (!conditional) { + _columnCategories.insert(column->GetName()); + } + } + auto &components = elem["components"].set_seq(); + // Components with renamed observables replace the original ones (Config::renameCombinedObservables). + auto renamed = _renamedComponents.find(simPdf.GetName()); + std::set replaced; + std::vector replacements; + // In the order of the columns and of their states, which is the order in which the categories are written to + // the domain. This does not depend on the order of the inputs of a RooSuperCategory, so the document is the same + // after a round trip. + for (std::vector const &key : selectionKeys(columns)) { + const std::string label = selectionLabel(indexCat, columns, key); + RooAbsPdf const *pdf = simPdf.getPdf(label); + if (!pdf) { + continue; + } + if (renamed != _renamedComponents.end()) { + auto found = renamed->second.find(label); + if (found != renamed->second.end()) { + replaced.insert(pdf); + pdf = found->second; + replacements.push_back(pdf); + } + } + JSONNode &component = components.append_child().set_map(); + component["key"].fill_seq(key); + component["dist"] << pdf->GetName(); + } + + std::vector servers; + for (RooAbsArg const *server : simPdf.servers()) { + if (replaced.count(server) == 0) { + servers.push_back(server); + } + } + servers.insert(servers.end(), replacements.begin(), replacements.end()); + exportObjects(servers, exportedObjectNames); +} + /** * @brief Import a function from the JSONNode into the workspace. * @@ -1354,6 +1711,23 @@ 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(); + } +} + +// A categorical axis: the labels in the order of the states of the category, which is the order of the bins of a +// RooDataHist. +void RooJSONFactoryWSTool::exportCategoricalAxis(JSONNode &obsNode, RooCategory const &cat) +{ + testValidName(cat.GetName(), false); + obsNode["name"] << cat.GetName(); + auto &labels = obsNode["categories"].set_seq(); + for (unsigned int i = 0; i < cat.size(); ++i) { + labels.append_child() << cat.getOrdinal(i).first; + } } /** @@ -1372,8 +1746,12 @@ void RooJSONFactoryWSTool::exportHisto(RooArgSet const &vars, std::size_t n, dou { auto &observablesNode = output["axes"].set_seq(); // axes have to be ordered to get consistent bin indices - for (auto *var : static_range_cast(vars)) { - exportAxis(observablesNode.append_child().set_map(), *var); + for (RooAbsArg const *var : vars) { + if (auto const *cat = dynamic_cast(var)) { + exportCategoricalAxis(observablesNode.append_child().set_map(), *cat); + } else { + exportAxis(observablesNode.append_child().set_map(), static_cast(*var)); + } } return exportArray(n, contents, output["contents"]); @@ -1420,83 +1798,234 @@ std::string makeValidNameOrError(std::string const &in) return out; } +// The labels of a selection key, made valid names and joined by "_". Used in the names of sub-data sets and of +// renamed observables. +std::string selectionKeyName(std::vector const &key) +{ + std::string out; + for (std::string const &label : key) { + out += (out.empty() ? "" : "_") + makeValidNameOrError(label); + } + return out; +} + +// The names of the category columns of `vars`. +std::vector selectionColumnNames(RooAbsCollection const &vars) +{ + std::vector names; + for (RooCategory const *column : selectionColumns(vars)) { + names.push_back(column->GetName()); + } + return names; +} + +// The RooSimultaneous of the workspace (the first by name) whose columns are `columns`, or nullptr. +RooSimultaneous const *findMatchingSimultaneous(RooWorkspace const &workspace, std::vector const &columns) +{ + std::vector candidates; + for (RooAbsArg *arg : workspace.allPdfs()) { + if (auto *candidate = dynamic_cast(arg)) { + candidates.push_back(candidate); + } + } + std::sort(candidates.begin(), candidates.end(), + [](auto *a, auto *b) { return std::string(a->GetName()) < std::string(b->GetName()); }); + for (RooSimultaneous const *candidate : candidates) { + std::vector candidateColumns; + for (RooCategory const *column : selectionColumns(candidate->indexCat())) { + candidateColumns.push_back(column->GetName()); + } + if (candidateColumns == columns) { + return candidate; + } + } + return nullptr; +} + } // namespace /** - * @brief Export a RooAbsCategory object to a JSONNode. + * @brief Export a data set with category observables as a "combined" data set. * - * This function exports a RooAbsCategory object, represented by the provided categories and indices, - * to a JSONNode. The category labels and corresponding indices are added to the JSONNode as key-value pairs. - * - * @param cat The RooAbsCategory object to be exported. - * @param node The JSONNode to which the category data will be exported. - * @return void + * The categories are the columns, sorted by name, and are declared as categorical axes of the default domain. There + * is one sub-data set per key, i.e. per combination of states, written like any other data set and named after the + * data set and the labels of the key. A key is written also without entries: RooFit treats a channel without entries + * as present and empty. If the workspace has a RooSimultaneous with the same columns, each sub-data set only carries + * the observables of the corresponding component, and keys without entries for which it has no component are left + * out; otherwise, every sub-data set carries all observables. */ -void RooJSONFactoryWSTool::exportCategory(RooAbsCategory const &cat, JSONNode &node) -{ - auto &labels = node["labels"].set_seq(); - auto &indices = node["indices"].set_seq(); - - for (auto const &item : cat) { - labels.append_child() << makeValidNameOrError(item.first); - indices.append_child() << item.second; +// Config::renameCombinedObservables: work out, before anything is exported, the new names of the observables of +// every combined data set, and the components of the matching RooSimultaneous with these observables. One +// observable per key is named "obs_", several are named "obs__". Observables that already +// have their new name are left alone, so that a second export of an imported document gives the same result. +void RooJSONFactoryWSTool::prepareCombinedObservableRenaming() +{ + _combinedObservableRenames.clear(); + _renamedComponents.clear(); + _renamedObjects.clear(); + if (!config().renameCombinedObservables) { + return; + } + std::vector datasets; + for (RooAbsData *data : _workspace.allData()) { + datasets.push_back(data); + } + sortByName(datasets); + for (RooAbsData const *data : datasets) { + std::vector columns = selectionColumns(*data->get()); + if (columns.empty()) { + continue; + } + RooSimultaneous const *simPdf = findMatchingSimultaneous(_workspace, selectionColumnNames(*data->get())); + for (std::vector const &key : selectionKeys(columns)) { + const std::string keyName = selectionKeyName(key); + RooArgSet candidates; + RooAbsPdf const *component = nullptr; + std::string label; + if (simPdf) { + label = selectionLabel(simPdf->indexCat(), columns, key); + component = simPdf->getPdf(label); + if (!component) { + continue; + } + std::unique_ptr observables{component->getObservables(*data->get())}; + candidates.add(*observables); + } else { + candidates.add(*data->get()); + } + std::vector observables; + for (RooAbsArg *arg : candidates) { + if (!dynamic_cast(arg)) { + observables.push_back(arg); + } + } + sortByName(observables); + std::map renames; + for (RooAbsArg const *obs : observables) { + std::string const newName = + "obs_" + keyName + (observables.size() == 1 ? "" : "_" + std::string(obs->GetName())); + if (newName == obs->GetName()) { + continue; + } + if (_workspace.arg(newName.c_str())) { + error("can't rename the observable '" + std::string(obs->GetName()) + "' of the data set '" + + data->GetName() + "' to '" + newName + "': an object of that name already exists"); + } + renames[obs->GetName()] = newName; + } + if (renames.empty()) { + continue; + } + _combinedObservableRenames[data->GetName()][keyName] = renames; + if (!simPdf || _renamedComponents[simPdf->GetName()].count(label) > 0) { + continue; + } + // The component with the renamed observables. Nodes that depend on them are cloned with the suffix + // "_", so that nodes shared by several components get one copy per component. + RooCustomizer customizer{*component, keyName.c_str()}; + for (auto const &[from, to] : renames) { + RooAbsArg *original = candidates.find(from.c_str()); + _renamedObjects.emplace_back(static_cast(original->clone(to.c_str()))); + customizer.replaceArg(*original, *_renamedObjects.back()); + } + auto *renamedComponent = static_cast(customizer.build()); + _renamedObjects.emplace_back(renamedComponent); + _renamedComponents[simPdf->GetName()][label] = renamedComponent; + } } } -// Split `data` by its index category into per-channel datasets and export each, returning the resulting -// component-name map. RooJSONFactoryWSTool::CombinedData RooJSONFactoryWSTool::exportCombinedData(RooAbsData const &data) { - RooAbsCategory *cat = findCategoryObservable(data); - - // prepare return value - RooJSONFactoryWSTool::CombinedData datamap; - - if (!cat) - return datamap; - // this is a combined dataset - - datamap.name = data.GetName(); - - // Write information necessary to reconstruct the combined dataset upon import - auto &child = getRooFitInternal(*_rootnodeOutput, "combined_datasets").set_map()[data.GetName()].set_map(); - child["index_cat"] << cat->GetName(); - exportCategory(*cat, child); - - // Find a RooSimultaneous model that would fit to this dataset - RooSimultaneous const *simPdf = nullptr; - auto *combinedPdfInfoNode = findRooFitInternal(*_rootnodeOutput, "combined_distributions"); - if (combinedPdfInfoNode) { - for (auto &info : combinedPdfInfoNode->children()) { - if (info["index_cat"].val() == cat->GetName()) { - simPdf = static_cast(_workspace.pdf(info.key())); + std::vector columns = selectionColumns(*data.get()); + CombinedData info; + info.name = data.GetName(); + for (RooCategory const *column : columns) { + info.columns.push_back(column->GetName()); + _domains->readCategory(*column); + } + + // Split the data by the key: by the index category of a matching RooSimultaneous, or by the columns themselves. + RooSimultaneous const *simPdf = findMatchingSimultaneous(_workspace, info.columns); + std::unique_ptr superCat; + RooAbsCategory const *indexCat = nullptr; + std::vector> parts; + if (simPdf) { + indexCat = &simPdf->indexCat(); + parts = data.split(*simPdf, true); + } else if (columns.size() == 1) { + indexCat = columns.front(); + parts = data.split(*indexCat, true); + } else { + RooArgSet inputs; + for (RooCategory const *column : columns) { + inputs.add(*column); + } + const std::string superName = info.name + "_key"; + superCat = std::make_unique(superName.c_str(), superName.c_str(), inputs); + indexCat = superCat.get(); + parts = data.split(*indexCat, true); + } + + // Fill the "combined" entry before exporting the sub-data sets: appending to "data" can invalidate references. + std::vector subDataSets; + { + JSONNode &output = appendNamedChild((*_rootnodeOutput)["data"], info.name); + output["type"] << "combined"; + output["columns"].fill_seq(info.columns); + auto &components = output["components"].set_seq(); + for (std::vector const &key : selectionKeys(columns)) { + std::string const label = selectionLabel(*indexCat, columns, key); + auto found = std::find_if(parts.begin(), parts.end(), [&](auto &d) { return d->GetName() == label; }); + if (found == parts.end() || ((*found)->numEntries() == 0 && simPdf && !simPdf->getPdf(label))) { + continue; } + const std::string keyName = selectionKeyName(key); + const std::string subName = info.name + "_" + keyName; + (*found)->SetName(subName.c_str()); + // Config::renameCombinedObservables: give the observables of this key their own names. + auto renamesOfData = _combinedObservableRenames.find(info.name); + if (renamesOfData != _combinedObservableRenames.end()) { + auto renames = renamesOfData->second.find(keyName); + if (renames != renamesOfData->second.end()) { + for (auto const &[from, to] : renames->second) { + (*found)->changeObservableName(from.c_str(), to.c_str()); + } + } + } + JSONNode &component = components.append_child().set_map(); + component["key"].fill_seq(key); + component["dataset"] << subName; + subDataSets.push_back(found->get()); } } - - // If there is an associated simultaneous pdf for the index category, we - // use the RooAbsData::split() overload that takes the RooSimultaneous. - // Like this, the observables that are not relevant for a given channel - // are automatically split from the component datasets. - std::vector> dataList{simPdf ? data.split(*simPdf, true) : data.split(*cat, true)}; - - for (std::unique_ptr const &absData : dataList) { - std::string catName(absData->GetName()); - std::string dataName = makeValidNameOrError(catName); - absData->SetName((std::string(data.GetName()) + "_" + dataName).c_str()); - datamap.components[catName] = absData->GetName(); - this->exportData(*absData); + // In HS3, the components of a combined distribution share no variate: a quantity measured in several components is a + // separate variate in each of them, e.g. "mass_SR" and "mass_CR". The sub-data sets, which are split by the + // components of the RooSimultaneous, must therefore not share observables. + if (simPdf) { + std::vector shared = sharedVariableNames(subDataSets, [](RooAbsData const *d) { return d; }); + if (!shared.empty()) { + error("the data set '" + info.name + "' is split by the components of '" + simPdf->GetName() + + "', which share the observables " + joinQuoted(shared) + + ". HS3 requires each component of a combined distribution to have its own variates: give the " + "observable a separate name in each component (e.g. 'mass_SR' and 'mass_CR' instead of 'mass')."); + } + } + for (RooAbsData *subData : subDataSets) { + exportData(*subData); } - return datamap; + return info; } // Export a single dataset `data` (binned or unbinned) to the output JSON. void RooJSONFactoryWSTool::exportData(RooAbsData const &data) { - RooAbsCategory *cat = findCategoryObservable(data); - - if (cat) - return; + for (RooAbsArg const *obs : *data.get()) { + if (dynamic_cast(obs) && !dynamic_cast(obs)) { + error("the data set '" + std::string(data.GetName()) + "' has the categorical observable '" + obs->GetName() + + "', which is not a RooCategory, and can't be written as a combined data set"); + } + } JSONNode &output = appendNamedChild((*_rootnodeOutput)["data"], data.GetName()); @@ -1512,10 +2041,19 @@ void RooJSONFactoryWSTool::exportData(RooAbsData const &data) } // this is a regular binned dataset + // Categories are categorical axes, declared by a categorical axis of the default domain. + auto readIntoDomain = [&](RooAbsArg const &var) { + if (auto const *cat = dynamic_cast(&var)) { + _domains->readCategory(*cat); + } else { + _domains->readVariable(static_cast(var)); + } + }; + if (auto dh = dynamic_cast(&data)) { output["type"] << "binned"; - for (auto *var : static_range_cast(variables)) { - _domains->readVariable(*var); + for (RooAbsArg const *var : variables) { + readIntoDomain(*var); } return exportHisto(variables, dh->numEntries(), dh->weightArray(), output); } @@ -1525,8 +2063,10 @@ void RooJSONFactoryWSTool::exportData(RooAbsData const &data) // RooDataSets from binned data to fit HistFactory models. In this case, it // doesn't make sense to export them like an unbinned dataset, because the // coordinates are redundant information with the binning. We only do this - // for 1D data for now. - if (data.isWeighted() && variables.size() == 1) { + // for 1D data for now. A variable without bins can't be binned data: unbinned data sets leave their axes without + // binning, and would otherwise be written as binned once they have no entries. + if (data.isWeighted() && variables.size() == 1 && dynamic_cast(variables[0]) && + static_cast(*variables[0]).getBins() > 0) { bool isBinnedData = false; auto &x = static_cast(*variables[0]); std::vector contents; @@ -1551,8 +2091,13 @@ void RooJSONFactoryWSTool::exportData(RooAbsData const &data) // this really is an unbinned dataset output["type"] << "unbinned"; auto &observablesNode = output["axes"].set_seq(); - for (auto *var : static_range_cast(variables)) { - _domains->readVariable(*var); + for (RooAbsArg const *arg : variables) { + readIntoDomain(*arg); + if (auto const *cat = dynamic_cast(arg)) { + exportCategoricalAxis(observablesNode.append_child().set_map(), *cat); + continue; + } + auto const *var = static_cast(arg); // HS3 forbids binning information on the axes of unbinned data. auto &axis = observablesNode.append_child().set_map(); testValidName(var->GetName(), false); @@ -1565,7 +2110,15 @@ void RooJSONFactoryWSTool::exportData(RooAbsData const &data) bool hasNonUnityWeights = false; for (int i = 0; i < data.numEntries(); ++i) { data.get(i); - coords.append_child().fill_seq(variables, [](auto x) { return static_cast(x)->getVal(); }); + // The entries of a categorical axis are the labels of the category. + auto &row = coords.append_child().set_seq(); + for (RooAbsArg const *arg : variables) { + if (auto const *cat = dynamic_cast(arg)) { + row.append_child() << cat->getCurrentLabel(); + } else { + row.append_child() << static_cast(arg)->getVal(); + } + } if (data.isWeighted()) { weightVals.push_back(data.weight()); if (data.weight() != 1.) @@ -1591,7 +2144,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 +2186,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 +2302,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. @@ -1725,40 +2378,40 @@ void RooJSONFactoryWSTool::exportModelConfig(JSONNode &rootnode, RooFit::ModelCo const std::vector &combDataSets, const std::vector &singleDataSets) { + // The analyses are named after the ModelConfig, so that ModelConfigs sharing a pdf and a data set (like an S+B + // and B-only pair) get separate analyses. auto pdf = mc.GetPdf(); auto simpdf = dynamic_cast(pdf); - if (simpdf) { - for (std::size_t i = 0; i < std::max(combDataSets.size(), std::size_t(1)); ++i) { - const bool hasdata = i < combDataSets.size(); - if (hasdata && !matches(combDataSets.at(i), simpdf)) - continue; - - std::string analysisName(mc.GetName()); - if (hasdata) - analysisName += "_" + combDataSets[i].name; - - exportSingleModelConfig(rootnode, mc, analysisName, hasdata ? &combDataSets[i].components : nullptr); + int founddata = 0; + if (simpdf && !simpdf->getAttribute(categoryConditionalTag)) { + // One analysis per combined data set with the same columns. + std::vector columns; + for (RooCategory const *column : selectionColumns(simpdf->indexCat())) { + columns.push_back(column->GetName()); + } + for (auto const &data : combDataSets) { + if (data.columns == columns) { + exportSingleModelConfig(rootnode, mc, std::string(mc.GetName()) + "_" + data.name, &data.name); + ++founddata; + } } } else { RooArgSet observables(*mc.GetObservables()); - int founddata = 0; for (auto *data : singleDataSets) { if (observables.equals(*(data->get()))) { - std::map mapping; - mapping[pdf->GetName()] = data->GetName(); - exportSingleModelConfig(rootnode, mc, std::string(mc.GetName()) + "_" + data->GetName(), &mapping); + const std::string dataName = data->GetName(); + exportSingleModelConfig(rootnode, mc, std::string(mc.GetName()) + "_" + dataName, &dataName); ++founddata; } } - if (founddata == 0) { - exportSingleModelConfig(rootnode, mc, mc.GetName(), nullptr); - } + } + if (founddata == 0) { + exportSingleModelConfig(rootnode, mc, mc.GetName(), nullptr); } } void RooJSONFactoryWSTool::exportSingleModelConfig(JSONNode &rootnode, RooFit::ModelConfig const &mc, - std::string const &analysisName, - std::map const *dataComponents) + std::string const &analysisName, std::string const *dataName) { auto pdf = mc.GetPdf(); @@ -1774,22 +2427,10 @@ void RooJSONFactoryWSTool::exportSingleModelConfig(JSONNode &rootnode, RooFit::M nllNode["distributions"].set_seq(); nllNode["data"].set_seq(); - if (dataComponents) { - auto simPdf = dynamic_cast(pdf); - if (simPdf) { - for (auto const &item : simPdf->indexCat()) { - const auto &dataComp = dataComponents->find(item.first); - nllNode["distributions"].append_child() << simPdf->getPdf(item.first)->GetName(); - nllNode["data"].append_child() << dataComp->second; - } - } else { - for (auto it : *dataComponents) { - nllNode["distributions"].append_child() << it.first; - nllNode["data"].append_child() << it.second; - } - } + nllNode["distributions"].append_child() << pdf->GetName(); + if (dataName) { + nllNode["data"].append_child() << *dataName; } else { - nllNode["distributions"].append_child() << pdf->GetName(); nllNode["data"].append_child() << 0; } @@ -1850,6 +2491,8 @@ void RooJSONFactoryWSTool::exportAllObjects(JSONNode &n) _domains = std::make_unique(); _varsNode = &makeVariablesNode(n); _rootnodeOutput = &n; + _columnCategories.clear(); + prepareCombinedObservableRenaming(); // export all toplevel pdfs std::vector allpdfs; @@ -1887,19 +2530,20 @@ void RooJSONFactoryWSTool::exportAllObjects(JSONNode &n) alldata.push_back(d); } sortByName(alldata); - // first, take care of combined datasets + // Data sets with category observables are combined data sets. std::vector singleData; std::vector combData; for (auto &d : alldata) { - auto data = this->exportCombinedData(*d); - if (!data.components.empty()) - combData.push_back(data); - else + // Without a matching RooSimultaneous, the categories of a data set are ordinary observables, written as + // categorical axes, unless its observables are renamed per key (Config::renameCombinedObservables). + if (selectionColumns(*d->get()).empty() || + (!config().renameCombinedObservables && + !findMatchingSimultaneous(_workspace, selectionColumnNames(*d->get())))) { singleData.push_back(d); - } - // next, take care datasets - for (auto &d : alldata) { - this->exportData(*d); + this->exportData(*d); + } else { + combData.push_back(this->exportCombinedData(*d)); + } } // export all ModelConfig objects and attached Pdfs @@ -2177,22 +2821,28 @@ void RooJSONFactoryWSTool::importAllNodes(const JSONNode &n) _attributesNode = nullptr; - // We delay the import of the data to after combineDatasets(), because it - // might be that some datasets are merged to combined datasets there. In - // that case, we will remove the components from the "datasets" vector so they - // don't get imported. + // The components of combined data sets are moved out of the "datasets" vector when the combined data set is + // assembled, so they don't get imported on their own. The ROOT-internal form of older files (combineDatasets()) + // is handled after the analyses, because their likelihoods refer to the components. std::vector> datasets; + std::map combinedData; if (auto dataNode = n.find("data")) { for (const auto &p : dataNode->children()) { + // Combined data sets are assembled from the others once all of them are loaded. + if (p.is_map() && p.has_child("type") && p["type"].val() == "combined") { + continue; + } datasets.push_back(loadData(p, _workspace, _domains.get())); } + combinedData = assembleCombinedData(*dataNode, _workspace, datasets); } // Now, read in analyses and likelihoods if there are any if (auto analysesNode = n.find("analyses")) { for (JSONNode const &analysisNode : analysesNode->children()) { - importAnalysis(*_rootnodeInput, analysisNode, n["likelihoods"], n["domains"], _workspace, datasets); + importAnalysis(*_rootnodeInput, analysisNode, n["likelihoods"], n["domains"], _workspace, datasets, + combinedData); } } diff --git a/roofit/hs3/test/testRooFitHS3.cxx b/roofit/hs3/test/testRooFitHS3.cxx index 560aea148a51b..d09f7b99c602f 100644 --- a/roofit/hs3/test/testRooFitHS3.cxx +++ b/roofit/hs3/test/testRooFitHS3.cxx @@ -9,9 +9,11 @@ #include #include #include +#include #include #include #include +#include #include #include #include @@ -35,6 +37,7 @@ #include #include #include +#include #include #include #include @@ -514,12 +517,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 +531,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 +882,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 +2639,8 @@ 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; + EXPECT_NE(exported.find(R"("constraints":["muConstraint"])"), 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 +2705,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":["alpha_lumiConstraint"])"), 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 +3348,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 +3358,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 +3405,22 @@ TEST(RooFitHS3, HistFactoryStatErrorConstraint) ASSERT_TRUE(RooJSONFactoryWSTool{ws}.importJSONfromString(makeJSON("Gauss", "Poisson"))); EXPECT_TRUE(hasGaussianStatConstraints(ws)); + // The constraints are exported as distributions 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(), "gaussian_dist") << json; + } RooWorkspace ws2; ASSERT_TRUE(RooJSONFactoryWSTool{ws2}.importJSONfromString(json)); @@ -3443,12 +3470,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 +3540,289 @@ 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 references the same constraints by name, and the result is stable under a round trip. + const std::string exported = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + EXPECT_NE(exported.find(R"("constraints":["alpha_a_constraint"])"), std::string::npos) << exported; + EXPECT_NE(exported.find(R"("constraints":["gamma_0_constraint",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 +} + +// 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 +3841,150 @@ 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":["g_c_0_constraint",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 @@ -3658,17 +4109,18 @@ TEST(RooFitHS3, ModelConfigsSharingPdfAndData) } { RooWorkspace ws; - ws.factory("Gaussian::pdf_a(x[0, -5, 5], mu[0, -5, 5], sigma[1, 0.1, 5])"); - ws.factory("Gaussian::pdf_b(x, mu, sigma)"); + // Each component has its own observable, as HS3 requires. + ws.factory("Gaussian::pdf_a(x_a[0, -5, 5], mu[0, -5, 5], sigma[1, 0.1, 5])"); + ws.factory("Gaussian::pdf_b(x_b[0, -5, 5], mu, sigma)"); ws.factory("SIMUL::sim_pdf(channel[a=0, b=1], a=pdf_a, b=pdf_b)"); - RooArgSet observables{*ws.var("x"), *ws.cat("channel")}; + RooArgSet observables{*ws.var("x_a"), *ws.var("x_b"), *ws.cat("channel")}; RooDataSet data{"data", "data", observables}; for (int index : {0, 1}) { ws.cat("channel")->setIndex(index); data.add(observables); } ws.import(data, RooFit::Silence()); - roundTrip(ws, "sim_pdf", "x,channel"); + roundTrip(ws, "sim_pdf", "x_a,x_b,channel"); } } @@ -3783,7 +4235,7 @@ 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":["alpha_sysConstraint"])"), std::string::npos) << json; RooWorkspace ws2; ASSERT_TRUE(RooJSONFactoryWSTool{ws2}.importJSONfromString(json)); @@ -4029,3 +4481,1202 @@ 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"); +} + +// A RooSimultaneous is a "combined_dist" (HS3 #119): the components are selected by a key with one label per column, +// and the columns are declared as categorical axes of the default domain. The ROOT-internal form of older files is +// still read. +TEST(RooFitHS3, CombinedDist) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + RooWorkspace ws; + ws.factory("Gaussian::model(x[-8, 8], mean[0, -8, 8], sigma[0.3, 0.1, 10])"); + ws.factory("Gaussian::model_ctl(x, mean, 1.0)"); + // The states are defined in this order, with non-contiguous indices. + ws.factory("SIMUL::simPdf(sample[physics=3, control=7], physics=model, control=model_ctl)"); + + const std::string json = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + auto tree = RooFit::Detail::JSONTree::create(json); + auto const &root = tree->rootnode(); + + // Structure of the written document. + auto const *simPdfNode = RooJSONFactoryWSTool::findNamedChild(root["distributions"], "simPdf"); + ASSERT_NE(simPdfNode, nullptr) << json; + EXPECT_EQ((*simPdfNode)["type"].val(), "combined_dist"); + ASSERT_EQ((*simPdfNode)["columns"].num_children(), 1u); + EXPECT_EQ((*simPdfNode)["columns"].child(0).val(), "sample"); + auto const &components = (*simPdfNode)["components"]; + ASSERT_EQ(components.num_children(), 2u); + EXPECT_EQ(components.child(0)["key"].num_children(), 1u); + EXPECT_EQ(components.child(0)["key"].child(0).val(), "physics"); + EXPECT_EQ(components.child(0)["dist"].val(), "model"); + EXPECT_EQ(components.child(1)["key"].child(0).val(), "control"); + EXPECT_EQ(components.child(1)["dist"].val(), "model_ctl"); + EXPECT_FALSE(root.find("misc", "ROOT_internal", "combined_distributions")) << json; + // The column is a categorical axis of the default domain, with the states in their order. + auto const *domain = RooJSONFactoryWSTool::findNamedChild(root["domains"], "default_domain"); + ASSERT_NE(domain, nullptr); + auto const *axis = RooJSONFactoryWSTool::findNamedChild((*domain)["axes"], "sample"); + ASSERT_NE(axis, nullptr) << json; + ASSERT_EQ((*axis)["categories"].num_children(), 2u); + EXPECT_EQ((*axis)["categories"].child(0).val(), "physics"); + EXPECT_EQ((*axis)["categories"].child(1).val(), "control"); + + // Round trip: the same RooSimultaneous, with the index category from the domain (indices = positions). + EXPECT_EQ(validate(ws, "simPdf"), 0); + RooWorkspace ws2; + ASSERT_TRUE(RooJSONFactoryWSTool{ws2}.importJSONfromString(json)); + auto *simPdf = dynamic_cast(ws2.pdf("simPdf")); + ASSERT_NE(simPdf, nullptr); + EXPECT_STREQ(simPdf->indexCat().GetName(), "sample"); + EXPECT_EQ(simPdf->indexCat().lookupIndex("physics"), 0); + EXPECT_EQ(simPdf->indexCat().lookupIndex("control"), 1); + EXPECT_STREQ(simPdf->getPdf("physics")->GetName(), "model"); + EXPECT_STREQ(simPdf->getPdf("control")->GetName(), "model_ctl"); + + // Files written by earlier ROOT versions describe the RooSimultaneous in misc.ROOT_internal. + RooWorkspace wsLegacy; + ASSERT_TRUE(RooJSONFactoryWSTool{wsLegacy}.importJSONfromString(R"({ + "metadata": {"hs3_version": "0.2"}, + "domains": [{"name": "default_domain", "type": "product_domain", + "axes": [{"name": "x", "min": -8.0, "max": 8.0}, {"name": "mean", "min": -8.0, "max": 8.0}]}], + "distributions": [ + {"name": "model", "type": "gaussian_dist", "x": "x", "mean": "mean", "sigma": 0.3}, + {"name": "model_ctl", "type": "gaussian_dist", "x": "x", "mean": "mean", "sigma": 1.0} + ], + "misc": {"ROOT_internal": {"combined_distributions": {"simPdf": { + "index_cat": "sample", "labels": ["physics", "control"], "indices": [3, 7], + "distributions": ["model", "model_ctl"]}}}} + })")); + auto *legacySimPdf = dynamic_cast(wsLegacy.pdf("simPdf")); + ASSERT_NE(legacySimPdf, nullptr); + EXPECT_EQ(legacySimPdf->indexCat().lookupIndex("control"), 7); + EXPECT_STREQ(legacySimPdf->getPdf("control")->GetName(), "model_ctl"); + // Its export is in the standard form. + const std::string fromLegacy = RooJSONFactoryWSTool{wsLegacy}.exportJSONtoString(); + EXPECT_NE(fromLegacy.find(R"("type":"combined_dist")"), std::string::npos) << fromLegacy; + EXPECT_EQ(fromLegacy.find("combined_distributions"), std::string::npos) << fromLegacy; +} + +// A RooSimultaneous over a RooSuperCategory has several columns, sorted by name, and keys with one label per column. +TEST(RooFitHS3, CombinedDistSeveralColumns) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + RooRealVar x{"x", "x", 0.0, 10.0}; + x.setBins(5); + RooGaussian sr{"sr", "sr", x, 3.0, 1.0}; + RooGaussian cr{"cr", "cr", x, 6.0, 1.0}; + RooGaussian srGlobal{"sr_global", "sr_global", x, 5.0, 2.0}; + // The super category is built over (region, kind); the columns are written sorted by name: (kind, region). + RooCategory region{"region", "region", {{"SR", 0}, {"CR", 1}}}; + RooCategory kind{"kind", "kind", {{"physics", 0}, {"global", 1}}}; + RooSuperCategory index{"index", "index", RooArgSet{region, kind}}; + // Three of the four keys have a component. + RooSimultaneous simPdf{ + "simPdf", "simPdf", {{"{SR;physics}", &sr}, {"{CR;physics}", &cr}, {"{SR;global}", &srGlobal}}, index}; + + RooWorkspace ws; + ws.import(simPdf, RooFit::Silence()); + const std::string json = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + auto tree = RooFit::Detail::JSONTree::create(json); + auto const *simPdfNode = RooJSONFactoryWSTool::findNamedChild(tree->rootnode()["distributions"], "simPdf"); + ASSERT_NE(simPdfNode, nullptr) << json; + ASSERT_EQ((*simPdfNode)["columns"].num_children(), 2u); + EXPECT_EQ((*simPdfNode)["columns"].child(0).val(), "kind"); + EXPECT_EQ((*simPdfNode)["columns"].child(1).val(), "region"); + std::map> keys; + for (auto const &component : (*simPdfNode)["components"].children()) { + ASSERT_EQ(component["key"].num_children(), 2u); + keys[component["dist"].val()] = {component["key"].child(0).val(), component["key"].child(1).val()}; + } + const std::map> expectedKeys{ + {"sr", {"physics", "SR"}}, {"cr", {"physics", "CR"}}, {"sr_global", {"global", "SR"}}}; + EXPECT_EQ(keys, expectedKeys) << json; + + // Round trip: the super category is rebuilt over the columns, and every component evaluates to the same values. + RooWorkspace ws2; + ASSERT_TRUE(RooJSONFactoryWSTool{ws2}.importJSONfromString(json)); + auto *simPdf2 = dynamic_cast(ws2.pdf("simPdf")); + ASSERT_NE(simPdf2, nullptr); + ASSERT_NE(dynamic_cast(&simPdf2->indexCat()), nullptr); + EXPECT_EQ(RooJSONFactoryWSTool{ws2}.exportJSONtoString(), json); + RooArgSet nset{x}; + RooArgSet nset2{*ws2.var("x")}; + for (auto const &item : expectedKeys) { + kind.setLabel(item.second[0]); + region.setLabel(item.second[1]); + ws2.cat("kind")->setLabel(item.second[0]); + ws2.cat("region")->setLabel(item.second[1]); + EXPECT_STREQ(simPdf2->getPdf(simPdf2->indexCat().getCurrentLabel())->GetName(), item.first.c_str()); + for (int i = 0; i < x.numBins(); ++i) { + x.setBin(i); + ws2.var("x")->setBin(i); + EXPECT_DOUBLE_EQ(simPdf2->getVal(nset2), simPdf.getVal(nset)) << item.first << " bin " << i; + } + } +} + +TEST(RooFitHS3, CombinedDistImportErrors) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + auto makeJSON = [](std::string const &columns, std::string const &components, + std::string const &axes = R"({"name": "sample", "categories": ["a", "b"]})") { + return R"({ + "metadata": {"hs3_version": "0.3"}, + "domains": [{"name": "default_domain", "type": "product_domain", + "axes": [{"name": "x", "min": 0.0, "max": 10.0}, )" + + axes + R"(]}], + "distributions": [ + {"name": "ga", "type": "gaussian_dist", "x": "x", "mean": 3.0, "sigma": 1.0}, + {"name": "gb", "type": "gaussian_dist", "x": "x", "mean": 6.0, "sigma": 1.0}, + {"name": "simPdf", "type": "combined_dist", "columns": )" + + columns + R"(, "components": )" + components + R"(} + ] + })"; + }; + const std::string good = R"([{"key": ["a"], "dist": "ga"}, {"key": ["b"], "dist": "gb"}])"; + { + RooWorkspace ws; + ASSERT_TRUE(RooJSONFactoryWSTool{ws}.importJSONfromString(makeJSON(R"(["sample"])", good))); + EXPECT_NE(dynamic_cast(ws.pdf("simPdf")), nullptr); + } + expectImportThrowsWithError( + makeJSON(R"(["channel"])", good), + "the column 'channel' of the combined_dist 'simPdf' is not declared by a categorical axis"); + expectImportThrowsWithError(makeJSON(R"([])", good), "must be a non-empty list of names of categorical axes"); + expectImportThrowsWithError(makeJSON(R"(["sample", "sample"])", R"([{"key": ["a", "a"], "dist": "ga"}])"), + "the column 'sample' of the combined_dist 'simPdf' is listed twice"); + expectImportThrowsWithError(makeJSON(R"(["sample"])", R"([{"key": ["a", "b"], "dist": "ga"}])"), + "must have a \"key\" with 1 labels"); + expectImportThrowsWithError(makeJSON(R"(["sample"])", R"([{"key": ["c"], "dist": "ga"}])"), + "the label 'c' in position 0 of a key of the combined_dist 'simPdf' is not one of the " + "categories of 'sample'"); + expectImportThrowsWithError( + makeJSON(R"(["sample"])", R"([{"key": ["a"], "dist": "ga"}, {"key": ["a"], "dist": "gb"}])"), + "the key [\"a\"] is used by two components"); + expectImportThrowsWithError(makeJSON(R"(["sample"])", R"([])"), + "\"components\" of the combined_dist 'simPdf' must be a non-empty list"); + expectImportThrowsWithError(makeJSON(R"(["sample"])", R"([{"key": ["a"]}])"), "must have a \"key\" and a \"dist\""); + // Whether a combined distribution is extended follows from its components: declaring it is not allowed, whatever + // the value. + expectImportThrowsWithError(makeJSON(R"(["sample"])", good + R"(, "extended": false)"), + "the combined_dist 'simPdf' must not have \"extended\""); + // The column must be categorical. + expectImportThrowsWithError( + makeJSON(R"(["sample"])", good, R"({"name": "sample", "min": 0.0, "max": 1.0})"), + "the column 'sample' of the combined_dist 'simPdf' is not declared by a categorical axis"); +} + +// A data set with a category observable is a "combined" data set (HS3 #119): one sub-data set per state, selected +// by a key, with the category declared as a categorical axis of the default domain. A likelihood pairs the combined +// data set with the combined_dist. The ROOT-internal form of older files is still read. +TEST(RooFitHS3, CombinedData) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + RooWorkspace ws{"ws"}; + ws.factory("Gaussian::model_1(x_1[0, 10], mean[5, 0, 10], sigma_1[0.5, 0.1, 10])"); + ws.factory("Gaussian::model_2(x_2[0, 10], mean, sigma_2[0.8, 0.1, 10])"); + ws.factory("SIMUL::simPdf(channelCat[channel_1=0, channel_2=1], channel_1=model_1, channel_2=model_2)"); + std::unique_ptr d1{ws.pdf("model_1")->generate(*ws.var("x_1"), 20)}; + std::unique_ptr d2{ws.pdf("model_2")->generate(*ws.var("x_2"), 30)}; + std::map parts{{"channel_1", d1.get()}, {"channel_2", d2.get()}}; + RooDataSet obsData{"obsData", + "obsData", + {*ws.var("x_1"), *ws.var("x_2")}, + RooFit::Index(*ws.cat("channelCat")), + RooFit::Import(parts)}; + ws.import(obsData); + RooFit::ModelConfig mc{"ModelConfig"}; + mc.SetWS(ws); + mc.SetPdf("simPdf"); + mc.SetParametersOfInterest("mean"); + mc.SetObservables("x_1,x_2,channelCat"); + ws.import(mc); + + const std::string json = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + auto tree = RooFit::Detail::JSONTree::create(json); + auto const &root = tree->rootnode(); + + // Structure: the combined entry, one sub-data set per state with only the observables of its component, the + // domain axis, the likelihood pairing the combined objects. + auto const *combined = RooJSONFactoryWSTool::findNamedChild(root["data"], "obsData"); + ASSERT_NE(combined, nullptr) << json; + EXPECT_EQ((*combined)["type"].val(), "combined"); + ASSERT_EQ((*combined)["columns"].num_children(), 1u); + EXPECT_EQ((*combined)["columns"].child(0).val(), "channelCat"); + EXPECT_FALSE(combined->has_child("axes")); + ASSERT_EQ((*combined)["components"].num_children(), 2u); + EXPECT_EQ((*combined)["components"].child(0)["key"].child(0).val(), "channel_1"); + EXPECT_EQ((*combined)["components"].child(0)["dataset"].val(), "obsData_channel_1"); + EXPECT_EQ((*combined)["components"].child(1)["key"].child(0).val(), "channel_2"); + EXPECT_EQ((*combined)["components"].child(1)["dataset"].val(), "obsData_channel_2"); + for (const char *channel : {"1", "2"}) { + auto const *part = RooJSONFactoryWSTool::findNamedChild(root["data"], std::string("obsData_channel_") + channel); + ASSERT_NE(part, nullptr) << json; + EXPECT_EQ((*part)["type"].val(), "unbinned"); + ASSERT_EQ((*part)["axes"].num_children(), 1u); + EXPECT_EQ((*part)["axes"].child(0)["name"].val(), std::string("x_") + channel); + } + EXPECT_FALSE(root.find("misc", "ROOT_internal", "combined_datasets")) << json; + auto const *domain = RooJSONFactoryWSTool::findNamedChild(root["domains"], "default_domain"); + ASSERT_NE(domain, nullptr); + EXPECT_NE(RooJSONFactoryWSTool::findNamedChild((*domain)["axes"], "channelCat"), nullptr) << json; + auto const *likelihood = RooJSONFactoryWSTool::findNamedChild(root["likelihoods"], "ModelConfig_obsData_likelihood"); + ASSERT_NE(likelihood, nullptr) << json; + ASSERT_EQ((*likelihood)["distributions"].num_children(), 1u); + EXPECT_EQ((*likelihood)["distributions"].child(0).val(), "simPdf"); + ASSERT_EQ((*likelihood)["data"].num_children(), 1u); + EXPECT_EQ((*likelihood)["data"].child(0).val(), "obsData"); + EXPECT_NE(RooJSONFactoryWSTool::findNamedChild(root["analyses"], "ModelConfig_obsData"), nullptr) << json; + + // Round trip: one data set with the category column, the same entries per state, the ModelConfig observes the + // category, and the document is stable. + RooWorkspace ws2; + ASSERT_TRUE(RooJSONFactoryWSTool{ws2}.importJSONfromString(json)); + EXPECT_EQ(ws2.allData().size(), 1u); + RooAbsData *obsData2 = ws2.data("obsData"); + ASSERT_NE(obsData2, nullptr); + EXPECT_EQ(obsData2->numEntries(), 50); + auto *channelCat2 = dynamic_cast(obsData2->get()->find("channelCat")); + ASSERT_NE(channelCat2, nullptr); + auto split2 = obsData2->split(*channelCat2); + ASSERT_EQ(split2.size(), 2u); + for (auto const &part : split2) { + EXPECT_EQ(part->numEntries(), std::string(part->GetName()) == "channel_1" ? 20 : 30) << part->GetName(); + } + auto *mc2 = dynamic_cast(ws2.obj("ModelConfig")); + ASSERT_NE(mc2, nullptr); + EXPECT_EQ(mc2->GetPdf(), ws2.pdf("simPdf")); + EXPECT_NE(mc2->GetObservables()->find("channelCat"), nullptr); + EXPECT_EQ(RooJSONFactoryWSTool{ws2}.exportJSONtoString(), json); + + // Without a RooSimultaneous over the category, the category is an ordinary observable: the data set is written + // with a categorical axis (HS3 combined data sets require a combined distribution to pair with), and read back + // with the same entries per state. + { + RooWorkspace wsDataOnly; + wsDataOnly.import(obsData); + const std::string dataOnly = RooJSONFactoryWSTool{wsDataOnly}.exportJSONtoString(); + auto dataOnlyTree = RooFit::Detail::JSONTree::create(dataOnly); + auto const &dataNode = dataOnlyTree->rootnode()["data"]; + EXPECT_EQ(RooJSONFactoryWSTool::findNamedChild(dataNode, "obsData_channel_1"), nullptr) << dataOnly; + auto const *single = RooJSONFactoryWSTool::findNamedChild(dataNode, "obsData"); + ASSERT_NE(single, nullptr) << dataOnly; + EXPECT_EQ((*single)["type"].val(), "unbinned") << dataOnly; + auto const *channelAxis = RooJSONFactoryWSTool::findNamedChild((*single)["axes"], "channelCat"); + ASSERT_NE(channelAxis, nullptr) << dataOnly; + EXPECT_EQ((*channelAxis)["categories"].num_children(), 2u) << dataOnly; + + RooWorkspace wsDataOnly2; + ASSERT_TRUE(RooJSONFactoryWSTool{wsDataOnly2}.importJSONfromString(dataOnly)); + RooAbsData *read = wsDataOnly2.data("obsData"); + ASSERT_NE(read, nullptr); + EXPECT_EQ(read->numEntries(), 50); + auto split = read->split(*static_cast(read->get()->find("channelCat"))); + ASSERT_EQ(split.size(), 2u); + for (auto const &part : split) { + EXPECT_EQ(part->numEntries(), std::string(part->GetName()) == "channel_1" ? 20 : 30) << part->GetName(); + } + EXPECT_EQ(RooJSONFactoryWSTool{wsDataOnly2}.exportJSONtoString(), dataOnly); + } + + // Files written by earlier ROOT versions: components named after the data set and the state, the combination + // described in misc.ROOT_internal, and the likelihood listing the components. + RooWorkspace wsLegacy; + ASSERT_TRUE(RooJSONFactoryWSTool{wsLegacy}.importJSONfromString(R"({ + "metadata": {"hs3_version": "0.2"}, + "domains": [{"name": "default_domain", "type": "product_domain", "axes": [ + {"name": "x_1", "min": 0.0, "max": 10.0}, {"name": "x_2", "min": 0.0, "max": 10.0}, {"name": "mean", "min": 0.0, "max": 10.0}]}], + "distributions": [ + {"name": "model_1", "type": "gaussian_dist", "x": "x_1", "mean": "mean", "sigma": 0.5}, + {"name": "model_2", "type": "gaussian_dist", "x": "x_2", "mean": "mean", "sigma": 0.8}], + "data": [ + {"name": "obsData_channel_1", "type": "unbinned", "axes": [{"name": "x_1", "min": 0.0, "max": 10.0}], "entries": [[1.0], [2.0]]}, + {"name": "obsData_channel_2", "type": "unbinned", "axes": [{"name": "x_2", "min": 0.0, "max": 10.0}], "entries": [[3.0]]}], + "likelihoods": [{"name": "simPdf_obsData_likelihood", "distributions": ["model_1", "model_2"], + "data": ["obsData_channel_1", "obsData_channel_2"]}], + "analyses": [{"name": "simPdf_obsData", "likelihood": "simPdf_obsData_likelihood", "parameters_of_interest": ["mean"]}], + "misc": {"ROOT_internal": { + "combined_datasets": {"obsData": {"index_cat": "channelCat", "labels": ["channel_1", "channel_2"], "indices": [0, 1]}}, + "combined_distributions": {"simPdf": {"index_cat": "channelCat", "labels": ["channel_1", "channel_2"], + "indices": [0, 1], "distributions": ["model_1", "model_2"]}}, + "ModelConfigs": {"simPdf_obsData": {"mcName": "ModelConfig", "pdfName": "simPdf"}}}} + })")); + EXPECT_EQ(wsLegacy.allData().size(), 1u); + ASSERT_NE(wsLegacy.data("obsData"), nullptr); + EXPECT_EQ(wsLegacy.data("obsData")->numEntries(), 3); + auto *mcLegacy = dynamic_cast(wsLegacy.obj("ModelConfig")); + ASSERT_NE(mcLegacy, nullptr); + EXPECT_EQ(mcLegacy->GetPdf(), wsLegacy.pdf("simPdf")); + const std::string fromLegacy = RooJSONFactoryWSTool{wsLegacy}.exportJSONtoString(); + EXPECT_NE(fromLegacy.find(R"("type":"combined")"), std::string::npos) << fromLegacy; + EXPECT_EQ(fromLegacy.find("combined_datasets"), std::string::npos) << fromLegacy; +} + +// A data set with several category observables has several columns, sorted by name. Keys without entries are +// written if the matching RooSimultaneous has a component for them, and left out otherwise. +TEST(RooFitHS3, CombinedDataSeveralColumns) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + // Each component has its own observable, as HS3 requires. + RooRealVar xSR{"x_sr", "x_sr", 0.0, 10.0}; + RooRealVar xCR{"x_cr", "x_cr", 0.0, 10.0}; + RooRealVar xGlobal{"x_sr_global", "x_sr_global", 0.0, 10.0}; + RooCategory region{"region", "region", {{"SR", 0}, {"CR", 1}}}; + RooCategory kind{"kind", "kind", {{"physics", 0}, {"global", 1}}}; + RooGaussian sr{"sr", "sr", xSR, 3.0, 1.0}; + RooGaussian cr{"cr", "cr", xCR, 6.0, 1.0}; + RooGaussian srGlobal{"sr_global", "sr_global", xGlobal, 5.0, 2.0}; + RooSuperCategory index{"index", "index", RooArgSet{region, kind}}; + RooSimultaneous simPdf{ + "simPdf", "simPdf", {{"{SR;physics}", &sr}, {"{CR;physics}", &cr}, {"{SR;global}", &srGlobal}}, index}; + // Entries per key; none for (physics, CR), and (global, CR) has no component. + const std::map, int> entries{{{"physics", "SR"}, 3}, {{"global", "SR"}, 1}}; + RooDataSet obs{"obs", "obs", {xSR, xCR, xGlobal, region, kind}}; + for (auto const &item : entries) { + kind.setLabel(item.first[0].c_str()); + region.setLabel(item.first[1].c_str()); + RooRealVar &x = item.first[0] == "global" ? xGlobal : xSR; + for (int i = 0; i < item.second; ++i) { + x.setVal(1.0 + i); + obs.add({xSR, xCR, xGlobal, region, kind}); + } + } + + RooWorkspace ws; + ws.import(simPdf, RooFit::Silence()); + ws.import(obs); + const std::string json = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + auto tree = RooFit::Detail::JSONTree::create(json); + auto const *combined = RooJSONFactoryWSTool::findNamedChild(tree->rootnode()["data"], "obs"); + ASSERT_NE(combined, nullptr) << json; + ASSERT_EQ((*combined)["columns"].num_children(), 2u); + EXPECT_EQ((*combined)["columns"].child(0).val(), "kind"); + EXPECT_EQ((*combined)["columns"].child(1).val(), "region"); + std::map, std::string> components; + for (auto const &component : (*combined)["components"].children()) { + components[{component["key"].child(0).val(), component["key"].child(1).val()}] = component["dataset"].val(); + } + const std::map, std::string> expected{{{"physics", "SR"}, "obs_physics_SR"}, + {{"physics", "CR"}, "obs_physics_CR"}, + {{"global", "SR"}, "obs_global_SR"}}; + EXPECT_EQ(components, expected) << json; + auto const *emptyPart = RooJSONFactoryWSTool::findNamedChild(tree->rootnode()["data"], "obs_physics_CR"); + ASSERT_NE(emptyPart, nullptr) << json; + EXPECT_EQ((*emptyPart)["entries"].num_children(), 0u); + + // Round trip: one data set with both category columns and the same entries per key; stable document. + RooWorkspace ws2; + ASSERT_TRUE(RooJSONFactoryWSTool{ws2}.importJSONfromString(json)); + EXPECT_EQ(ws2.allData().size(), 1u); + RooAbsData *obs2 = ws2.data("obs"); + ASSERT_NE(obs2, nullptr); + EXPECT_EQ(RooJSONFactoryWSTool{ws2}.exportJSONtoString(), json); + std::map, int> entries2; + for (int i = 0; i < obs2->numEntries(); ++i) { + RooArgSet const &row = *obs2->get(i); + ++entries2[{static_cast(row["kind"]).getCurrentLabel(), + static_cast(row["region"]).getCurrentLabel()}]; + } + EXPECT_EQ(entries2, entries); + // The RooSimultaneous splits it by its RooSuperCategory, which is rebuilt over the columns. + auto parts = obs2->split(*dynamic_cast(ws2.pdf("simPdf")), /*createEmptyDataSets=*/true); + std::map partEntries; + for (auto const &part : parts) { + partEntries[part->GetName()] = part->numEntries(); + } + const std::map expectedParts{ + {"{physics;SR}", 3}, {"{physics;CR}", 0}, {"{global;SR}", 1}, {"{global;CR}", 0}}; + EXPECT_EQ(partEntries, expectedParts); +} + +TEST(RooFitHS3, CombinedDataImportErrors) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + auto makeJSON = [](std::string const &combined, std::string const &axes = "") { + return R"({ + "metadata": {"hs3_version": "0.3"}, + "domains": [{"name": "default_domain", "type": "product_domain", + "axes": [{"name": "x", "min": 0.0, "max": 10.0}, {"name": "sample", "categories": ["a", "b"]}]}], + "data": [ + {"name": "da", "type": "unbinned", "axes": [{"name": "x", "min": 0.0, "max": 10.0}], "entries": [[1.0]]}, + {"name": "db", "type": "unbinned", "axes": [{"name": "x", "min": 0.0, "max": 10.0}], "entries": [[2.0], [3.0]]}, + {"name": "obs", "type": "combined", )" + + axes + combined + R"(} + ] + })"; + }; + const std::string good = + R"("columns": ["sample"], "components": [{"key": ["a"], "dataset": "da"}, {"key": ["b"], "dataset": "db"}])"; + { + RooWorkspace ws; + ASSERT_TRUE(RooJSONFactoryWSTool{ws}.importJSONfromString(makeJSON(good))); + EXPECT_EQ(ws.allData().size(), 1u); + ASSERT_NE(ws.data("obs"), nullptr); + EXPECT_EQ(ws.data("obs")->numEntries(), 3); + } + expectImportThrowsWithError(makeJSON(good, R"("axes": [{"name": "x", "min": 0.0, "max": 10.0}], )"), + "the combined data set 'obs' must not have \"axes\""); + expectImportThrowsWithError( + makeJSON(R"("columns": ["channel"], "components": [{"key": ["a"], "dataset": "da"}])"), + "the column 'channel' of the combined data set 'obs' is not declared by a categorical axis"); + expectImportThrowsWithError(makeJSON(R"("columns": ["sample"], "components": [])"), + "\"components\" of the combined data set 'obs' must be a non-empty list"); + expectImportThrowsWithError(makeJSON(R"("columns": ["sample"], "components": [{"key": ["a"]}])"), + "must have a \"key\" and a \"dataset\""); + expectImportThrowsWithError(makeJSON(R"("columns": ["sample"], "components": [{"key": ["c"], "dataset": "da"}])"), + "the label 'c' in position 0 of a key of the combined data set 'obs'"); + expectImportThrowsWithError(makeJSON(R"("columns": ["sample"], "components": [{"key": ["a"], "dataset": "dc"}])"), + "the component 'dc' of the combined data set 'obs' is not a data set of the document"); + expectImportThrowsWithError( + makeJSON( + R"("columns": ["sample"], "components": [{"key": ["a"], "dataset": "da"}, {"key": ["b"], "dataset": "da"}])"), + "the component 'da' of the combined data set 'obs' is not a data set of the document, or it is already a " + "component"); + expectImportThrowsWithError( + makeJSON( + R"("columns": ["sample"], "components": [{"key": ["a"], "dataset": "da"}, {"key": ["a"], "dataset": "db"}])"), + "the key [\"a\"] is used by two components of the combined data set 'obs'"); +} + +// A likelihood that pairs a combined_dist with a combined data set: HS3 treats a key without sub-data set as an +// unobserved component, which RooFit cannot represent (it would be a channel without entries, adding the expected +// events of an extended component to the likelihood). This is an error for extended components and a warning +// otherwise; a sub-data set without component is ignored with a warning. +TEST(RooFitHS3, CombinedDataPairing) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + // `components` of the distribution and of the data, and whether the components are extended. + auto makeJSON = [](std::string const &pdfComponents, std::string const &dataComponents, bool extended) { + auto component = [&](std::string const &label, double mean) { + const std::string gauss = R"({"name": "g_)" + label + R"(", "type": "gaussian_dist", "x": "x_)" + label + + R"(", "mean": )" + std::to_string(mean) + R"(, "sigma": 1.0})"; + if (!extended) { + return gauss; + } + return gauss + R"(, {"name": "e_)" + label + R"(", "type": "rate_extended_dist", "dist": "g_)" + label + + R"(", "rate": "n_)" + label + R"("})"; + }; + auto data = [](std::string const &label) { + return R"({"name": "d_)" + label + R"(", "type": "unbinned", "axes": [{"name": "x_)" + label + + R"(", "min": 0.0, "max": 10.0}], "entries": [[5.0]]})"; + }; + return R"({ + "metadata": {"hs3_version": "0.3"}, + "domains": [{"name": "default_domain", "type": "product_domain", "axes": [ + {"name": "x_a", "min": 0.0, "max": 10.0}, {"name": "x_b", "min": 0.0, "max": 10.0}, + {"name": "x_c", "min": 0.0, "max": 10.0}, {"name": "n_a", "min": 0.0, "max": 100.0}, + {"name": "n_b", "min": 0.0, "max": 100.0}, {"name": "n_c", "min": 0.0, "max": 100.0}, + {"name": "sample", "categories": ["a", "b", "c"]}]}], + "distributions": [)" + + component("a", 3.0) + ", " + component("b", 4.0) + ", " + component("c", 5.0) + R"(, + {"name": "simPdf", "type": "combined_dist", "columns": ["sample"], "components": )" + + pdfComponents + R"(}], + "data": [)" + + data("a") + ", " + data("b") + ", " + data("c") + R"(, + {"name": "obs", "type": "combined", "columns": ["sample"], "components": )" + + dataComponents + R"(}], + "likelihoods": [{"name": "nll", "distributions": ["simPdf"], "data": ["obs"]}], + "analyses": [{"name": "fit", "likelihood": "nll", "domains": ["default_domain"]}] + })"; + }; + auto pdfComponents = [](std::string const &prefix, std::vector const &labels) { + std::string out = "["; + for (auto const &label : labels) { + out += std::string(out.size() > 1 ? ", " : "") + R"({"key": [")" + label + R"("], "dist": ")" + prefix + + label + "\"}"; + } + return out + "]"; + }; + auto dataComponents = [](std::vector const &labels) { + std::string out = "["; + for (auto const &label : labels) { + out += + std::string(out.size() > 1 ? ", " : "") + R"({"key": [")" + label + R"("], "dataset": "d_)" + label + "\"}"; + } + return out + "]"; + }; + + // Every component has a sub-data set: no messages. + { + RooWorkspace ws; + const std::string warnings = captureMessages(RooFit::WARNING, RooFit::ObjectHandling, [&] { + ASSERT_TRUE(RooJSONFactoryWSTool{ws}.importJSONfromString( + makeJSON(pdfComponents("e_", {"a", "b", "c"}), dataComponents({"a", "b", "c"}), true))); + }); + EXPECT_EQ(warnings, ""); + auto *mc = dynamic_cast(ws.obj("fit")); + ASSERT_NE(mc, nullptr); + EXPECT_EQ(mc->GetPdf(), ws.pdf("simPdf")); + } + // An extended component without sub-data set is an error. + expectImportThrowsWithError(makeJSON(pdfComponents("e_", {"a", "b", "c"}), dataComponents({"a", "b"}), true), + "the combined data set 'obs' has no sub-data set for the key [\"c\"] of the extended " + "component 'e_c' of 'simPdf'"); + // A non-extended component without sub-data set is a warning, and so is a sub-data set without component. + { + RooWorkspace ws; + const std::string warnings = captureMessages(RooFit::WARNING, RooFit::ObjectHandling, [&] { + ASSERT_TRUE(RooJSONFactoryWSTool{ws}.importJSONfromString( + makeJSON(pdfComponents("g_", {"a", "b", "c"}), dataComponents({"a", "b"}), false))); + }); + EXPECT_NE(warnings.find("the combined data set 'obs' has no sub-data set for the key [\"c\"] of the component " + "'g_c' of 'simPdf', which RooFit treats as present and empty"), + std::string::npos) + << warnings; + } + { + RooWorkspace ws; + const std::string warnings = captureMessages(RooFit::WARNING, RooFit::ObjectHandling, [&] { + ASSERT_TRUE(RooJSONFactoryWSTool{ws}.importJSONfromString( + makeJSON(pdfComponents("e_", {"a", "b"}), dataComponents({"a", "b", "c"}), true))); + }); + EXPECT_NE(warnings.find("the sub-data set for the key [\"c\"] of the combined data set 'obs' is ignored, " + "because 'simPdf' has no component for it"), + std::string::npos) + << warnings; + } +} + +// A "category_conditional_dist" (HS3 #119) selects a component by the value of categorical parameters. It is a +// RooSimultaneous whose index category is a parameter: RooFit then evaluates the selected component. The +// RooSimultaneous carries a tag so that it is written back as a category_conditional_dist, and the value of the +// parameter is written to the parameter points. +TEST(RooFitHS3, CategoryConditionalDist) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + const std::string json = R"({ + "metadata": {"hs3_version": "0.3"}, + "domains": [{"name": "default_domain", "type": "product_domain", "axes": [ + {"name": "x", "min": 0.0, "max": 10.0}, {"name": "bkg_function", "categories": ["narrow", "wide"]}]}], + "parameter_points": [{"name": "point", "parameters": [{"name": "bkg_function", "value": "wide"}]}], + "distributions": [ + {"name": "narrow", "type": "gaussian_dist", "x": "x", "mean": 5.0, "sigma": 1.0}, + {"name": "wide", "type": "gaussian_dist", "x": "x", "mean": 5.0, "sigma": 3.0}, + {"name": "bkg", "type": "category_conditional_dist", "categorical_parameter": ["bkg_function"], + "components": [{"key": ["narrow"], "dist": "narrow"}, {"key": ["wide"], "dist": "wide"}]}] + })"; + RooWorkspace ws; + ASSERT_TRUE(RooJSONFactoryWSTool{ws}.importJSONfromString(json)); + auto *bkg = dynamic_cast(ws.pdf("bkg")); + ASSERT_NE(bkg, nullptr); + EXPECT_TRUE(bkg->getAttribute("hs3_category_conditional_dist")); + RooCategory *choice = ws.cat("bkg_function"); + ASSERT_NE(choice, nullptr); + EXPECT_EQ(&bkg->indexCat(), choice); + + // The value is the one of the selected component, conditional on the parameter (not in the normalization set). + RooRealVar &x = *ws.var("x"); + RooArgSet nset{x}; + ASSERT_TRUE(ws.loadSnapshot("point")); + EXPECT_STREQ(choice->getCurrentLabel(), "wide"); + for (double value : {2.0, 5.0, 8.0}) { + x.setVal(value); + choice->setLabel("wide"); + EXPECT_DOUBLE_EQ(bkg->getVal(nset), ws.pdf("wide")->getVal(nset)) << "x = " << value; + choice->setLabel("narrow"); + EXPECT_DOUBLE_EQ(bkg->getVal(nset), ws.pdf("narrow")->getVal(nset)) << "x = " << value; + } + + // Export: the same type, the parameter in the domain and in the parameter point, and a stable document. + ws.loadSnapshot("point"); + const std::string exported = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + auto tree = RooFit::Detail::JSONTree::create(exported); + auto const &root = tree->rootnode(); + auto const *bkgNode = RooJSONFactoryWSTool::findNamedChild(root["distributions"], "bkg"); + ASSERT_NE(bkgNode, nullptr) << exported; + EXPECT_EQ((*bkgNode)["type"].val(), "category_conditional_dist"); + ASSERT_EQ((*bkgNode)["categorical_parameter"].num_children(), 1u); + EXPECT_EQ((*bkgNode)["categorical_parameter"].child(0).val(), "bkg_function"); + EXPECT_FALSE(bkgNode->has_child("columns")); + auto const *point = RooJSONFactoryWSTool::findNamedChild(root["parameter_points"], "point"); + ASSERT_NE(point, nullptr) << exported; + auto const *choiceInPoint = RooJSONFactoryWSTool::findNamedChild((*point)["parameters"], "bkg_function"); + ASSERT_NE(choiceInPoint, nullptr) << exported; + EXPECT_EQ((*choiceInPoint)["value"].val(), "wide"); + auto const *domain = RooJSONFactoryWSTool::findNamedChild(root["domains"], "default_domain"); + ASSERT_NE(domain, nullptr); + EXPECT_NE(RooJSONFactoryWSTool::findNamedChild((*domain)["axes"], "bkg_function"), nullptr) << exported; + RooWorkspace ws2; + ASSERT_TRUE(RooJSONFactoryWSTool{ws2}.importJSONfromString(exported)); + ws2.loadSnapshot("point"); + EXPECT_EQ(RooJSONFactoryWSTool{ws2}.exportJSONtoString(), exported); + + // A RooSimultaneous built in RooFit is a combined_dist, unless it carries the tag. + { + RooWorkspace wsRooFit; + wsRooFit.factory("Gaussian::narrow(x[0, 10], 5.0, 1.0)"); + wsRooFit.factory("Gaussian::wide(x, 5.0, 3.0)"); + wsRooFit.factory("SIMUL::bkg(bkg_function[narrow=0, wide=1], narrow=narrow, wide=wide)"); + EXPECT_NE(RooJSONFactoryWSTool{wsRooFit}.exportJSONtoString().find(R"("type":"combined_dist")"), + std::string::npos); + wsRooFit.pdf("bkg")->setAttribute("hs3_category_conditional_dist"); + EXPECT_NE(RooJSONFactoryWSTool{wsRooFit}.exportJSONtoString().find(R"("type":"category_conditional_dist")"), + std::string::npos); + } + + // The parameter must be declared by a categorical axis of a domain. + expectImportThrowsWithError(R"({ + "metadata": {"hs3_version": "0.3"}, + "domains": [{"name": "default_domain", "type": "product_domain", "axes": [{"name": "x", "min": 0.0, "max": 10.0}]}], + "distributions": [ + {"name": "narrow", "type": "gaussian_dist", "x": "x", "mean": 5.0, "sigma": 1.0}, + {"name": "bkg", "type": "category_conditional_dist", "categorical_parameter": ["bkg_function"], + "components": [{"key": ["narrow"], "dist": "narrow"}]}] + })", + "the column 'bkg_function' of the category_conditional_dist 'bkg' is not declared by a " + "categorical axis of a domain"); +} + +// The names of a RooSimultaneous and of a combined data set are part of the standard document now, not of +// misc.ROOT_internal: a document without "misc" imports with the same names, and the ModelConfig finds its pdf +// through the likelihood. +TEST(RooFitHS3, CombinedObjectsKeepNamesWithoutMisc) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + // Remove the top-level key `key` (and the comma before it) from a compact JSON document. + auto removeTopLevelKey = [](std::string json, std::string const &key) { + const std::string needle = ",\"" + key + "\":"; + const std::size_t begin = json.find(needle); + if (begin == std::string::npos) { + return json; + } + std::size_t pos = begin + needle.size(); + int depth = 0; + bool inString = false; + for (; pos < json.size(); ++pos) { + const char c = json[pos]; + if (inString) { + if (c == '\\') { + ++pos; + } else if (c == '"') { + inString = false; + } + } else if (c == '"') { + inString = true; + } else if (c == '{' || c == '[') { + ++depth; + } else if (c == '}' || c == ']') { + if (--depth == 0) { + ++pos; + break; + } + } + } + return json.erase(begin, pos - begin); + }; + + RooWorkspace ws{"ws"}; + ws.factory("Gaussian::model_1(x_1[0, 10], mean[5, 0, 10], sigma_1[0.5, 0.1, 10])"); + ws.factory("Gaussian::model_2(x_2[0, 10], mean, sigma_2[0.8, 0.1, 10])"); + // Not named "simPdf", which older documents relied on as the default name. + ws.factory("SIMUL::channels(region[SR=0, CR=1], SR=model_1, CR=model_2)"); + std::unique_ptr d1{ws.pdf("model_1")->generate(*ws.var("x_1"), 20)}; + std::unique_ptr d2{ws.pdf("model_2")->generate(*ws.var("x_2"), 30)}; + std::map parts{{"SR", d1.get()}, {"CR", d2.get()}}; + RooDataSet observed{"observed", + "observed", + {*ws.var("x_1"), *ws.var("x_2")}, + RooFit::Index(*ws.cat("region")), + RooFit::Import(parts)}; + ws.import(observed); + RooFit::ModelConfig mc{"ModelConfig"}; + mc.SetWS(ws); + mc.SetPdf("channels"); + mc.SetParametersOfInterest("mean"); + mc.SetObservables("x_1,x_2,region"); + ws.import(mc); + + const std::string json = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + const std::string withoutMisc = removeTopLevelKey(json, "misc"); + ASSERT_NE(withoutMisc, json); + EXPECT_EQ(withoutMisc.find("\"misc\""), std::string::npos) << withoutMisc; + + RooWorkspace ws2; + ASSERT_TRUE(RooJSONFactoryWSTool{ws2}.importJSONfromString(withoutMisc)); + auto *channels = dynamic_cast(ws2.pdf("channels")); + ASSERT_NE(channels, nullptr); + EXPECT_STREQ(channels->indexCat().GetName(), "region"); + EXPECT_STREQ(channels->getPdf("SR")->GetName(), "model_1"); + EXPECT_STREQ(channels->getPdf("CR")->GetName(), "model_2"); + EXPECT_EQ(ws2.allData().size(), 1u); + RooAbsData *observed2 = ws2.data("observed"); + ASSERT_NE(observed2, nullptr); + EXPECT_EQ(observed2->numEntries(), 50); + EXPECT_NE(observed2->get()->find("region"), nullptr); + // Without misc.ROOT_internal.ModelConfigs, the ModelConfig is named after the analysis. + auto *mc2 = dynamic_cast(ws2.obj("ModelConfig_observed")); + ASSERT_NE(mc2, nullptr); + EXPECT_EQ(mc2->GetPdf(), channels); + EXPECT_NE(mc2->GetObservables()->find("region"), nullptr); + EXPECT_NE(mc2->GetParametersOfInterest()->find("mean"), nullptr); + + // The analysis is named after the ModelConfig, which got the name of the analysis: the names of the analysis and + // the likelihood change once on the re-export. From then on, the document is stable. + const std::string json2 = RooJSONFactoryWSTool{ws2}.exportJSONtoString(); + RooWorkspace ws3; + ASSERT_TRUE(RooJSONFactoryWSTool{ws3}.importJSONfromString(json2)); + EXPECT_EQ(RooJSONFactoryWSTool{ws3}.exportJSONtoString(), json2); +} + +// In HS3, the components of a combined distribution share no variate: a quantity measured in several components is a +// separate variate in each of them ("mass_SR", "mass_CR"). A RooSimultaneous whose components share an observable can +// therefore not be written with its data set. A document whose sub-data sets share an axis is still read, with a +// warning. +TEST(RooFitHS3, CombinedComponentsMustNotShareObservables) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + RooRealVar mass{"mass", "mass", 0.0, 10.0}; + RooCategory channel{"channel", "channel", {{"SR", 0}, {"CR", 1}}}; + RooGaussian sr{"sr", "sr", mass, 3.0, 1.0}; + RooGaussian cr{"cr", "cr", mass, 6.0, 1.0}; + RooSimultaneous simPdf{"simPdf", "simPdf", {{"SR", &sr}, {"CR", &cr}}, channel}; + RooDataSet obs{"obs", "obs", {mass, channel}}; + for (const char *label : {"SR", "CR"}) { + channel.setLabel(label); + mass.setVal(4.0); + obs.add({mass, channel}); + } + RooWorkspace ws; + ws.import(simPdf, RooFit::Silence()); + ws.import(obs); + const std::string errors = captureMessages(RooFit::ERROR, RooFit::IO, [&] { + EXPECT_THROW(RooJSONFactoryWSTool{ws}.exportJSONtoString(), std::runtime_error); + }); + EXPECT_NE(errors.find("which share the observables 'mass'"), std::string::npos) << errors; + EXPECT_NE(errors.find("'mass_SR' and 'mass_CR'"), std::string::npos) << errors; + + const std::string json = R"({ + "metadata": {"hs3_version": "0.3"}, + "domains": [{"name": "default_domain", "type": "product_domain", "axes": [ + {"name": "channel", "categories": ["SR", "CR"]}, {"name": "mass", "min": 0.0, "max": 10.0}]}], + "data": [ + {"name": "obs_SR", "type": "unbinned", "axes": [{"name": "mass", "min": 0.0, "max": 10.0}], "entries": [[4.0]]}, + {"name": "obs_CR", "type": "unbinned", "axes": [{"name": "mass", "min": 0.0, "max": 10.0}], "entries": [[5.0]]}, + {"name": "obs", "type": "combined", "columns": ["channel"], + "components": [{"key": ["SR"], "dataset": "obs_SR"}, {"key": ["CR"], "dataset": "obs_CR"}]} + ] + })"; + RooWorkspace imported; + const std::string warnings = captureMessages(RooFit::WARNING, RooFit::ObjectHandling, [&] { + EXPECT_TRUE(RooJSONFactoryWSTool{imported}.importJSONfromString(json)); + }); + EXPECT_NE(warnings.find("share the axes 'mass'"), std::string::npos) << warnings; + RooAbsData *data = imported.data("obs"); + ASSERT_NE(data, nullptr); + EXPECT_EQ(data->numEntries(), 2); +} + +// With RooJSONFactoryWSTool::config().renameCombinedObservables, the observables of combined data sets are renamed +// per key on export ("obs_SR", or "obs_SR_x" with several observables), and the components of the matching +// RooSimultaneous are written with them; nodes that depend on a renamed observable get one copy per component. This +// allows exporting a RooSimultaneous whose components share an observable. The likelihood stays the same, and a +// second export gives the same document. +TEST(RooFitHS3, RenameCombinedObservables) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + struct FlagGuard { + bool old = RooJSONFactoryWSTool::config().renameCombinedObservables; + FlagGuard() { RooJSONFactoryWSTool::config().renameCombinedObservables = true; } + ~FlagGuard() { RooJSONFactoryWSTool::config().renameCombinedObservables = old; } + } flagGuard; + + RooRealVar mass{"mass", "mass", 0.0, 10.0}; + RooCategory channel{"channel", "channel", {{"SR", 0}, {"CR", 1}}}; + // A function of the observable, shared by both components. + RooFormulaVar shift{"shift", "mass - 1.0", {mass}}; + RooGaussian sr{"sr", "sr", shift, 3.0, 1.0}; + RooGaussian cr{"cr", "cr", shift, 6.0, 1.5}; + RooSimultaneous simPdf{"simPdf", "simPdf", {{"SR", &sr}, {"CR", &cr}}, channel}; + RooDataSet obs{"obs", "obs", {mass, channel}}; + for (auto const &[label, value] : + std::vector>{{"SR", 3.5}, {"SR", 4.5}, {"CR", 6.0}, {"CR", 7.5}, {"CR", 8.0}}) { + channel.setLabel(label); + mass.setVal(value); + obs.add({mass, channel}); + } + std::unique_ptr nll{simPdf.createNLL(obs)}; + const double expectedNll = nll->getVal(); + + RooWorkspace ws; + ws.import(simPdf, RooFit::Silence()); + ws.import(obs); + const std::string json = RooJSONFactoryWSTool{ws}.exportJSONtoString(); + auto tree = RooFit::Detail::JSONTree::create(json); + auto const &root = tree->rootnode(); + + // The components and the functions of the observable, one per key, with the renamed observables. + auto const *combined = RooJSONFactoryWSTool::findNamedChild(root["distributions"], "simPdf"); + ASSERT_NE(combined, nullptr) << json; + std::map components; + for (auto const &component : (*combined)["components"].children()) { + components[component["key"].child(0).val()] = component["dist"].val(); + } + EXPECT_EQ(components, (std::map{{"SR", "sr_SR"}, {"CR", "cr_CR"}})) << json; + for (const char *key : {"SR", "CR"}) { + auto const *func = RooJSONFactoryWSTool::findNamedChild(root["functions"], std::string("shift_") + key); + ASSERT_NE(func, nullptr) << json; + EXPECT_NE((*func)["expression"].val().find(std::string("obs_") + key), std::string::npos) << json; + auto const *part = RooJSONFactoryWSTool::findNamedChild(root["data"], std::string("obs_") + key); + ASSERT_NE(part, nullptr) << json; + ASSERT_EQ((*part)["axes"].num_children(), 1u) << json; + EXPECT_EQ((*part)["axes"].child(0)["name"].val(), std::string("obs_") + key) << json; + } + EXPECT_EQ(RooJSONFactoryWSTool::findNamedChild(root["functions"], "shift"), nullptr) << json; + + // Read back, the likelihood is the same. + RooWorkspace imported; + ASSERT_TRUE(RooJSONFactoryWSTool{imported}.importJSONfromString(json)); + auto *simPdf2 = dynamic_cast(imported.pdf("simPdf")); + ASSERT_NE(simPdf2, nullptr); + RooAbsData *obs2 = imported.data("obs"); + ASSERT_NE(obs2, nullptr); + std::unique_ptr nll2{simPdf2->createNLL(*obs2)}; + EXPECT_NEAR(nll2->getVal(), expectedNll, 1e-10); + + // The observables already have their new names, so a second export gives the same document. + EXPECT_EQ(RooJSONFactoryWSTool{imported}.exportJSONtoString(), json); + + // Several observables per key: "obs__". Without a RooSimultaneous, all observables of the data + // set are renamed. + RooRealVar x{"x", "x", 0.0, 1.0}; + RooRealVar y{"y", "y", 0.0, 1.0}; + RooDataSet xy{"xy", "xy", {x, y, channel}}; + for (const char *label : {"SR", "CR"}) { + channel.setLabel(label); + xy.add({x, y, channel}); + } + RooWorkspace wsXY; + wsXY.import(xy); + const std::string jsonXY = RooJSONFactoryWSTool{wsXY}.exportJSONtoString(); + auto treeXY = RooFit::Detail::JSONTree::create(jsonXY); + auto const *partSR = RooJSONFactoryWSTool::findNamedChild(treeXY->rootnode()["data"], "xy_SR"); + ASSERT_NE(partSR, nullptr) << jsonXY; + std::set axes; + for (auto const &axis : (*partSR)["axes"].children()) { + axes.insert(axis["name"].val()); + } + EXPECT_EQ(axes, (std::set{"obs_SR_x", "obs_SR_y"})) << jsonXY; +} 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" } }