Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
JSONFactories_HistFactory.cxx
Go to the documentation of this file.
1/*
2 * Project: RooFit
3 * Authors:
4 * Carsten D. Burgard, DESY/ATLAS, Dec 2021
5 *
6 * Copyright (c) 2022, CERN
7 *
8 * Redistribution and use in source and binary forms,
9 * with or without modification, are permitted according to the terms
10 * listed in LICENSE (http://roofit.sourceforge.net/license.txt)
11 */
12
14#include <RooFitHS3/JSONIO.h>
16
21#include <RooConstVar.h>
22#include <RooRealVar.h>
23#include <RooDataHist.h>
24#include <RooHistFunc.h>
25#include <RooRealSumPdf.h>
26#include <RooBinWidthFunction.h>
27#include <RooProdPdf.h>
28#include <RooPoisson.h>
29#include <RooFormulaVar.h>
30#include <RooLognormal.h>
31#include <RooGaussian.h>
32#include <RooProduct.h>
33#include <RooWorkspace.h>
34
35#include <regex>
36
37#include "static_execute.h"
38#include "JSONIOUtils.h"
39
41
42using namespace RooStats::HistFactory;
43using namespace RooStats::HistFactory::Detail;
45
46namespace {
47
48inline void writeAxis(JSONNode &axis, RooRealVar const &obs)
49{
50 auto &binning = obs.getBinning();
51 if (binning.isUniform()) {
52 axis["nbins"] << obs.numBins();
53 axis["min"] << obs.getMin();
54 axis["max"] << obs.getMax();
55 } else {
56 auto &edges = axis["edges"];
57 edges.set_seq();
58 double val = binning.binLow(0);
59 edges.append_child() << val;
60 for (int i = 0; i < binning.numBins(); ++i) {
61 val = binning.binHigh(i);
62 edges.append_child() << val;
63 }
64 }
65}
66
67double round_prec(double d, int nSig)
68{
69 if (d == 0.0)
70 return 0.0;
71 int ndigits = std::floor(std::log10(std::abs(d))) + 1 - nSig;
72 double sf = std::pow(10, ndigits);
73 if (std::abs(d / sf) < 2)
74 ndigits--;
75 return sf * std::round(d / sf);
76}
77
78// To avoid repeating the same string literals that can potentially get out of
79// sync.
80namespace Literals {
81constexpr auto staterror = "staterror";
82}
83
84void erasePrefix(std::string &str, std::string_view prefix)
85{
86 if (startsWith(str, prefix)) {
87 str.erase(0, prefix.size());
88 }
89}
90
91bool eraseSuffix(std::string &str, std::string_view suffix)
92{
93 if (endsWith(str, suffix)) {
94 str.erase(str.size() - suffix.size());
95 return true;
96 } else {
97 return false;
98 }
99}
100
101template <class Coll>
102void sortByName(Coll &coll)
103{
104 std::sort(coll.begin(), coll.end(), [](auto &l, auto &r) { return l.name < r.name; });
105}
106
107template <class T>
108T *findClient(RooAbsArg *gamma)
109{
110 for (const auto &client : gamma->clients()) {
111 if (auto casted = dynamic_cast<T *>(client)) {
112 return casted;
113 } else {
114 T *c = findClient<T>(client);
115 if (c)
116 return c;
117 }
118 }
119 return nullptr;
120}
121
123{
124 if (!g)
125 return nullptr;
127 if (constraint_p)
128 return constraint_p;
130 if (constraint_g)
131 return constraint_g;
133 if (constraint_l)
134 return constraint_l;
135 return nullptr;
136}
137
138std::string toString(TClass *c)
139{
140 if (!c) {
141 return "Const";
142 }
143 if (c == RooPoisson::Class()) {
144 return "Poisson";
145 }
146 if (c == RooGaussian::Class()) {
147 return "Gauss";
148 }
149 if (c == RooLognormal::Class()) {
150 return "Lognormal";
151 }
152 return "unknown";
153}
154
155inline std::string defaultGammaName(std::string const &sysname, std::size_t i)
156{
157 return "gamma_" + sysname + "_bin_" + std::to_string(i);
158}
159
160/// Export the names of the gamma parameters to the modifier struct if the
161/// names don't match the default gamma parameter names, which is gamma_<sysname>_bin_<i>
162void optionallyExportGammaParameters(JSONNode &mod, std::string const &sysname, std::vector<RooAbsReal *> const &params,
163 bool forceExport = true)
164{
165 std::vector<std::string> paramNames;
166 bool needExport = forceExport;
167 for (std::size_t i = 0; i < params.size(); ++i) {
168 std::string name(params[i]->GetName());
169 paramNames.push_back(name);
170 if (name != defaultGammaName(sysname, i)) {
171 needExport = true;
172 }
173 }
174 if (needExport) {
175 mod["parameters"].fill_seq(paramNames);
176 }
177}
178
179RooRealVar &createNominal(RooWorkspace &ws, std::string const &parname, double val, double min, double max)
180{
181 RooRealVar &nom = getOrCreate<RooRealVar>(ws, "nom_" + parname, val, min, max);
182 nom.setConstant(true);
183 return nom;
184}
185
186/// Get the conventional name of the constraint pdf for a constrained
187/// parameter.
188std::string constraintName(std::string const &paramName)
189{
190 return paramName + "Constraint";
191}
192
193ParamHistFunc &createPHF(const std::string &phfname, std::string const &sysname,
194 const std::vector<std::string> &parnames, const std::vector<double> &vals,
195 RooJSONFactoryWSTool &tool, RooAbsCollection &constraints, const RooArgSet &observables,
196 const std::string &constraintType, double gammaMin, double gammaMax, double minSigma)
197{
198 RooWorkspace &ws = *tool.workspace();
199
201 for (std::size_t i = 0; i < vals.size(); ++i) {
202 const std::string name = parnames.empty() ? defaultGammaName(sysname, i) : parnames[i];
204 }
205
206 auto &phf = tool.wsEmplace<ParamHistFunc>(phfname, observables, gammas);
207
208 if (constraintType != "Const") {
210 gammas, vals, minSigma, constraintType == "Poisson" ? Constraint::Poisson : Constraint::Gaussian);
211 for (auto const &term : constraintsInfo.constraints) {
213 constraints.add(*ws.pdf(term->GetName()));
214 }
215 } else {
216 for (auto *gamma : static_range_cast<RooRealVar *>(gammas)) {
217 gamma->setConstant(true);
218 }
219 }
220
221 return phf;
222}
223
224bool hasStaterror(const JSONNode &comp)
225{
226 if (!comp.has_child("modifiers"))
227 return false;
228 for (const auto &mod : comp["modifiers"].children()) {
229 if (mod["type"].val() == ::Literals::staterror)
230 return true;
231 }
232 return false;
233}
234
235const JSONNode &findStaterror(const JSONNode &comp)
236{
237 if (comp.has_child("modifiers")) {
238 for (const auto &mod : comp["modifiers"].children()) {
239 if (mod["type"].val() == ::Literals::staterror)
240 return mod;
241 }
242 }
243 RooJSONFactoryWSTool::error("sample '" + RooJSONFactoryWSTool::name(comp) + "' does not have a " +
244 ::Literals::staterror + " modifier!");
245}
246
247RooAbsPdf &
248getOrCreateConstraint(RooJSONFactoryWSTool &tool, const JSONNode &mod, RooRealVar &param, const std::string &sample)
249{
250 if (auto constrName = mod.find("constraint_name")) {
251 auto constraint_name = constrName->val();
252 auto constraint = tool.workspace()->pdf(constraint_name);
253 if (!constraint) {
254 constraint = tool.request<RooAbsPdf>(constrName->val(), sample);
255 }
256 if (!constraint) {
257 RooJSONFactoryWSTool::error("unable to find definition of of constraint '" + constraint_name +
258 "' for modifier '" + RooJSONFactoryWSTool::name(mod) + "'");
259 }
260 if (auto gauss = dynamic_cast<RooGaussian *const>(constraint)) {
261 param.setError(gauss->getSigma().getVal());
262 }
263 return *constraint;
264 } else {
265 std::cout << "creating new constraint for " << param << std::endl;
266 std::string constraint_type = "Gauss";
267 if (auto constrType = mod.find("constraint_type")) {
269 }
270 if (constraint_type == "Gauss") {
271 param.setError(1.0);
272 return getOrCreate<RooGaussian>(*tool.workspace(), constraintName(param.GetName()), param,
273 *tool.workspace()->var(std::string("nom_") + param.GetName()), 1.);
274 }
275 RooJSONFactoryWSTool::error("unknown or invalid constraint for modifier '" + RooJSONFactoryWSTool::name(mod) +
276 "'");
277 }
278}
279
281 RooAbsArg const *mcStatObject, const std::string &fprefix, const JSONNode &p,
282 RooArgSet &constraints)
283{
284 RooWorkspace &ws = *tool.workspace();
285
287 std::string prefixedName = fprefix + "_" + sampleName;
288
289 std::string channelName = fprefix;
290 erasePrefix(channelName, "model_");
291
292 if (!p.has_child("data")) {
293 RooJSONFactoryWSTool::error("sample '" + sampleName + "' does not define a 'data' key");
294 }
295
296 auto &hf = tool.wsEmplace<RooHistFunc>("hist_" + prefixedName, varlist, dh);
297 hf.SetTitle(RooJSONFactoryWSTool::name(p).c_str());
298
301
302 shapeElems.add(tool.wsEmplace<RooBinWidthFunction>(prefixedName + "_binWidth", hf, true));
303
304 if (hasStaterror(p)) {
306 }
307
308 if (p.has_child("modifiers")) {
310 std::vector<double> overall_low;
311 std::vector<double> overall_high;
312 std::vector<int> overall_interp;
313
317
318 int idx = 0;
319 for (const auto &mod : p["modifiers"].children()) {
320 std::string const &modtype = mod["type"].val();
321 std::string const &sysname =
322 mod.has_child("name")
323 ? mod["name"].val()
324 : (mod.has_child("parameter") ? mod["parameter"].val() : "syst_" + std::to_string(idx));
325 ++idx;
326 if (modtype == "staterror") {
327 // this is dealt with at a different place, ignore it for now
328 } else if (modtype == "normfactor") {
331 if (mod.has_child("constraint_name") || mod.has_child("constraint_type")) {
332 // for norm factors, constraints are optional
334 }
335 } else if (modtype == "normsys") {
336 auto *parameter = mod.find("parameter");
337 std::string parname(parameter ? parameter->val() : "alpha_" + sysname);
338 createNominal(ws, parname, 0.0, -10, 10);
339 auto &par = getOrCreate<RooRealVar>(ws, parname, 0., -5, 5);
340 overall_nps.add(par);
341 auto &data = mod["data"];
342 int interp = 4;
343 if (mod.has_child("interpolation")) {
344 interp = mod["interpolation"].val_int();
345 }
346 double low = data["lo"].val_double();
347 double high = data["hi"].val_double();
348
349 // the below contains a a hack to cut off variations that go below 0
350 // this is needed because with interpolation code 4, which is the default, interpolation is done in
351 // log-space. hence, values <= 0 result in NaN which propagate throughout the model and cause evaluations to
352 // fail if you know a nicer way to solve this, please go ahead and fix the lines below
353 if (interp == 4 && low <= 0)
354 low = std::numeric_limits<double>::epsilon();
355 if (interp == 4 && high <= 0)
356 high = std::numeric_limits<double>::epsilon();
357
358 overall_low.push_back(low);
359 overall_high.push_back(high);
360 overall_interp.push_back(interp);
361
362 constraints.add(getOrCreateConstraint(tool, mod, par, sampleName));
363 } else if (modtype == "histosys") {
364 auto *parameter = mod.find("parameter");
365 std::string parname(parameter ? parameter->val() : "alpha_" + sysname);
366 createNominal(ws, parname, 0.0, -10, 10);
367 auto &par = getOrCreate<RooRealVar>(ws, parname, 0., -5, 5);
368 histNps.add(par);
369 auto &data = mod["data"];
370 histoLo.add(tool.wsEmplace<RooHistFunc>(
371 sysname + "Low_" + prefixedName, varlist,
373 histoHi.add(tool.wsEmplace<RooHistFunc>(
374 sysname + "High_" + prefixedName, varlist,
375 RooJSONFactoryWSTool::readBinnedData(data["hi"], sysname + "High_" + prefixedName, varlist)));
376 constraints.add(getOrCreateConstraint(tool, mod, par, sampleName));
377 } else if (modtype == "shapesys") {
378 std::string funcName = channelName + "_" + sysname + "_ShapeSys";
379 // funcName should be "<channel_name>_<sysname>_ShapeSys"
380 std::vector<double> vals;
381 for (const auto &v : mod["data"]["vals"].children()) {
382 vals.push_back(v.val_double());
383 }
384 std::vector<std::string> parnames;
385 for (const auto &v : mod["parameters"].children()) {
386 parnames.push_back(v.val());
387 }
388 if (vals.empty()) {
389 RooJSONFactoryWSTool::error("unable to instantiate shapesys '" + sysname + "' with 0 values!");
390 }
391 std::string constraint(mod.has_child("constraint_type") ? mod["constraint_type"].val()
392 : mod.has_child("constraint") ? mod["constraint"].val()
393 : "unknown");
394 shapeElems.add(createPHF(funcName, sysname, parnames, vals, tool, constraints, varlist, constraint,
396 } else if (modtype == "custom") {
397 RooAbsReal *obj = ws.function(sysname);
398 if (!obj) {
399 RooJSONFactoryWSTool::error("unable to find custom modifier '" + sysname + "'");
400 }
401 if (obj->dependsOn(varlist)) {
402 shapeElems.add(*obj);
403 } else {
404 normElems.add(*obj);
405 }
406 } else {
407 RooJSONFactoryWSTool::error("modifier '" + sysname + "' of unknown type '" + modtype + "'");
408 }
409 }
410
411 std::string interpName = sampleName + "_" + channelName + "_epsilon";
412 if (!overall_nps.empty()) {
415 normElems.add(v);
416 }
417 if (!histNps.empty()) {
418 auto &v = tool.wsEmplace<PiecewiseInterpolation>("histoSys_" + prefixedName, hf, histoLo, histoHi, histNps);
420 v.setAllInterpCodes(4); // default interpCode for HistFactory
421 shapeElems.add(v);
422 } else {
423 shapeElems.add(hf);
424 }
425 }
426
427 tool.wsEmplace<RooProduct>(prefixedName + "_shapes", shapeElems);
428 if (!normElems.empty()) {
429 tool.wsEmplace<RooProduct>(prefixedName + "_scaleFactors", normElems);
430 } else {
431 ws.factory("RooConstVar::" + prefixedName + "_scaleFactors(1.)");
432 }
433
434 return true;
435}
436
437class HistFactoryImporter : public RooFit::JSONIO::Importer {
438public:
439 bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override
440 {
441 std::string name = RooJSONFactoryWSTool::name(p);
442 if (!p.has_child("samples")) {
443 RooJSONFactoryWSTool::error("no samples in '" + name + "', skipping.");
444 }
445 double statErrThresh = 0;
446 std::string statErrType = "Poisson";
447 if (p.has_child(::Literals::staterror)) {
448 auto &staterr = p[::Literals::staterror];
449 if (staterr.has_child("relThreshold"))
450 statErrThresh = staterr["relThreshold"].val_double();
451 if (staterr.has_child("constraint_type"))
452 statErrType = staterr["constraint_type"].val();
453 }
454 std::vector<double> sumW;
455 std::vector<double> sumW2;
456 std::vector<std::string> gammaParnames;
458
459 std::string fprefix = name;
460
461 std::vector<std::unique_ptr<RooDataHist>> data;
462 for (const auto &comp : p["samples"].children()) {
463 std::unique_ptr<RooDataHist> dh = RooJSONFactoryWSTool::readBinnedData(
464 comp["data"], fprefix + "_" + RooJSONFactoryWSTool::name(comp) + "_dataHist", observables);
465 size_t nbins = dh->numEntries();
466
467 if (hasStaterror(comp)) {
468 if (sumW.empty()) {
469 sumW.resize(nbins);
470 sumW2.resize(nbins);
471 }
472 for (size_t i = 0; i < nbins; ++i) {
473 sumW[i] += dh->weight(i);
474 sumW2[i] += dh->weightSquared(i);
475 }
476 if (gammaParnames.empty()) {
477 if (auto staterrorParams = findStaterror(comp).find("parameters")) {
478 for (const auto &v : staterrorParams->children()) {
479 gammaParnames.push_back(v.val());
480 }
481 }
482 }
483 }
484 data.emplace_back(std::move(dh));
485 }
486
487 RooAbsArg *mcStatObject = nullptr;
488 RooArgSet constraints;
489 if (!sumW.empty()) {
490 std::string channelName = name;
491 erasePrefix(channelName, "model_");
492
493 std::vector<double> errs(sumW.size());
494 for (size_t i = 0; i < sumW.size(); ++i) {
495 errs[i] = std::sqrt(sumW2[i]) / sumW[i];
496 // avoid negative sigma. This NP will be set constant anyway later
497 errs[i] = std::max(errs[i], 0.);
498 }
499
501 &createPHF("mc_stat_" + channelName, "stat_" + channelName, gammaParnames, errs, *tool, constraints,
503 }
504
505 int idx = 0;
507 RooArgList coefs;
508 for (const auto &comp : p["samples"].children()) {
509 importHistSample(*tool, *data[idx], observables, mcStatObject, fprefix, comp, constraints);
510 ++idx;
511
512 std::string const &compName = RooJSONFactoryWSTool::name(comp);
513 funcs.add(*tool->request<RooAbsReal>(fprefix + "_" + compName + "_shapes", name));
514 coefs.add(*tool->request<RooAbsReal>(fprefix + "_" + compName + "_scaleFactors", name));
515 }
516
517 if (constraints.empty()) {
518 tool->wsEmplace<RooRealSumPdf>(name, funcs, coefs, true);
519 } else {
520 std::string sumName = name + "_model";
521 erasePrefix(sumName, "model_");
522 auto &sum = tool->wsEmplace<RooRealSumPdf>(sumName, funcs, coefs, true);
523 sum.SetTitle(name.c_str());
524 tool->wsEmplace<RooProdPdf>(name, constraints, RooFit::Conditional(sum, observables));
525 }
526 return true;
527 }
528};
529
530class FlexibleInterpVarStreamer : public RooFit::JSONIO::Exporter {
531public:
532 std::string const &key() const override
533 {
534 static const std::string keystring = "interpolation0d";
535 return keystring;
536 }
537 bool exportObject(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &elem) const override
538 {
539 auto fip = static_cast<const RooStats::HistFactory::FlexibleInterpVar *>(func);
540 elem["type"] << key();
541 elem["interpolationCodes"].fill_seq(fip->interpolationCodes());
542 RooJSONFactoryWSTool::fillSeq(elem["vars"], fip->variables());
543 elem["nom"] << fip->nominal();
544 elem["high"].fill_seq(fip->high(), fip->variables().size());
545 elem["low"].fill_seq(fip->low(), fip->variables().size());
546 return true;
547 }
548};
549
550class PiecewiseInterpolationStreamer : public RooFit::JSONIO::Exporter {
551public:
552 std::string const &key() const override
553 {
554 static const std::string keystring = "interpolation";
555 return keystring;
556 }
557 bool exportObject(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &elem) const override
558 {
559 const PiecewiseInterpolation *pip = static_cast<const PiecewiseInterpolation *>(func);
560 elem["type"] << key();
561 elem["interpolationCodes"].fill_seq(pip->interpolationCodes());
562 elem["positiveDefinite"] << pip->positiveDefinite();
563 RooJSONFactoryWSTool::fillSeq(elem["vars"], pip->paramList());
564 elem["nom"] << pip->nominalHist()->GetName();
565 RooJSONFactoryWSTool::fillSeq(elem["high"], pip->highList(), pip->paramList().size());
566 RooJSONFactoryWSTool::fillSeq(elem["low"], pip->lowList(), pip->paramList().size());
567 return true;
568 }
569};
570
571class PiecewiseInterpolationFactory : public RooFit::JSONIO::Importer {
572public:
573 bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override
574 {
575 std::string name(RooJSONFactoryWSTool::name(p));
576
577 RooArgList vars{tool->requestArgList<RooRealVar>(p, "vars")};
578
579 auto &pip = tool->wsEmplace<PiecewiseInterpolation>(name, *tool->requestArg<RooAbsReal>(p, "nom"),
580 tool->requestArgList<RooAbsReal>(p, "low"),
581 tool->requestArgList<RooAbsReal>(p, "high"), vars);
582
583 pip.setPositiveDefinite(p["positiveDefinite"].val_bool());
584
585 if (p.has_child("interpolationCodes")) {
586 std::size_t i = 0;
587 for (auto const &node : p["interpolationCodes"].children()) {
588 pip.setInterpCode(*static_cast<RooAbsReal *>(vars.at(i)), node.val_int(), true);
589 ++i;
590 }
591 }
592
593 return true;
594 }
595};
596
597class FlexibleInterpVarFactory : public RooFit::JSONIO::Importer {
598public:
599 bool importArg(RooJSONFactoryWSTool *tool, const JSONNode &p) const override
600 {
601 std::string name(RooJSONFactoryWSTool::name(p));
602 if (!p.has_child("high")) {
603 RooJSONFactoryWSTool::error("no high variations of '" + name + "'");
604 }
605 if (!p.has_child("low")) {
606 RooJSONFactoryWSTool::error("no low variations of '" + name + "'");
607 }
608 if (!p.has_child("nom")) {
609 RooJSONFactoryWSTool::error("no nominal variation of '" + name + "'");
610 }
611
612 double nom(p["nom"].val_double());
613
614 RooArgList vars{tool->requestArgList<RooRealVar>(p, "vars")};
615
616 std::vector<double> high;
617 high << p["high"];
618
619 std::vector<double> low;
620 low << p["low"];
621
622 if (vars.size() != low.size() || vars.size() != high.size()) {
623 RooJSONFactoryWSTool::error("FlexibleInterpVar '" + name +
624 "' has non-matching lengths of 'vars', 'high' and 'low'!");
625 }
626
627 auto &fip = tool->wsEmplace<RooStats::HistFactory::FlexibleInterpVar>(name, vars, nom, low, high);
628
629 if (p.has_child("interpolationCodes")) {
630 size_t i = 0;
631 for (auto const &node : p["interpolationCodes"].children()) {
632 fip.setInterpCode(*static_cast<RooAbsReal *>(vars.at(i)), node.val_int());
633 ++i;
634 }
635 }
636
637 return true;
638 }
639};
640
641struct NormFactor {
642 std::string name;
643 RooAbsReal const *param = nullptr;
644 RooAbsPdf const *constraint = nullptr;
645 TClass *constraintType = RooGaussian::Class();
646 NormFactor(RooAbsReal const &par, const RooAbsPdf *constr = nullptr)
647 : name{par.GetName()}, param{&par}, constraint{constr}
648 {
649 }
650};
651
652struct NormSys {
653 std::string name = "";
654 RooAbsReal const *param = nullptr;
655 double low = 1.;
656 double high = 1.;
657 int interpolationCode = 4;
658 RooAbsPdf const *constraint = nullptr;
659 TClass *constraintType = RooGaussian::Class();
660 NormSys() {};
661 NormSys(const std::string &n, RooAbsReal *const p, double h, double l, int i, const RooAbsPdf *c)
662 : name(n), param(p), low(l), high(h), interpolationCode(i), constraint(c), constraintType(c->IsA())
663 {
664 }
665};
666
667struct HistoSys {
668 std::string name;
669 RooAbsReal const *param = nullptr;
670 std::vector<double> low;
671 std::vector<double> high;
672 RooAbsPdf const *constraint = nullptr;
673 TClass *constraintType = RooGaussian::Class();
674 HistoSys(const std::string &n, RooAbsReal *const p, RooHistFunc *l, RooHistFunc *h, const RooAbsPdf *c)
675 : name(n), param(p), constraint(c), constraintType(c->IsA())
676 {
677 low.assign(l->dataHist().weightArray(), l->dataHist().weightArray() + l->dataHist().numEntries());
678 high.assign(h->dataHist().weightArray(), h->dataHist().weightArray() + h->dataHist().numEntries());
679 }
680};
681struct ShapeSys {
682 std::string name;
683 std::vector<double> constraints;
684 std::vector<RooAbsReal *> parameters;
685 RooAbsPdf const *constraint = nullptr;
686 TClass *constraintType = RooGaussian::Class();
687 ShapeSys(const std::string &n) : name{n} {}
688};
689
690struct GenericElement {
691 std::string name;
692 RooAbsReal *function = nullptr;
693 GenericElement(RooAbsReal *e) : name(e->GetName()), function(e) {};
694};
695
696std::string stripOuterParens(const std::string &s)
697{
698 size_t start = 0;
699 size_t end = s.size();
700
701 while (start < end && s[start] == '(' && s[end - 1] == ')') {
702 int depth = 0;
703 bool balanced = true;
704 for (size_t i = start; i < end - 1; ++i) {
705 if (s[i] == '(')
706 ++depth;
707 else if (s[i] == ')')
708 --depth;
709 if (depth == 0 && i < end - 1) {
710 balanced = false;
711 break;
712 }
713 }
714 if (balanced) {
715 ++start;
716 --end;
717 } else {
718 break;
719 }
720 }
721 return s.substr(start, end - start);
722}
723
724std::vector<std::string> splitTopLevelProduct(const std::string &expr)
725{
726 std::vector<std::string> parts;
727 int depth = 0;
728 size_t start = 0;
729 bool foundTopLevelStar = false;
730
731 for (size_t i = 0; i < expr.size(); ++i) {
732 char c = expr[i];
733 if (c == '(') {
734 ++depth;
735 } else if (c == ')') {
736 --depth;
737 } else if (c == '*' && depth == 0) {
738 foundTopLevelStar = true;
739 std::string sub = expr.substr(start, i - start);
740 parts.push_back(stripOuterParens(sub));
741 start = i + 1;
742 }
743 }
744
745 if (!foundTopLevelStar) {
746 return {}; // Not a top-level product
747 }
748
749 std::string sub = expr.substr(start);
750 parts.push_back(stripOuterParens(sub));
751 return parts;
752}
753
754#include <regex>
755#include <string>
756#include <cctype>
757#include <cstdlib>
758#include <iostream>
759
760NormSys parseOverallModifierFormula(const std::string &s, RooFormulaVar *formula)
761{
762 static const std::regex pattern(
763 R"(^\s*1(?:\.0)?\s*([\+\-])\s*([a-zA-Z_][a-zA-Z0-9_]*|[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)\s*\*\s*([a-zA-Z_][a-zA-Z0-9_]*|[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)\s*$)");
764
765 NormSys sys;
766 double sign = 1.0;
767
768 std::smatch match;
769 if (std::regex_match(s, match, pattern)) {
770 if (match[1].str() == "-") {
771 sign = -1.0;
772 }
773
774 std::string token2 = match[2].str();
775 std::string token3 = match[4].str();
776
777 RooAbsReal *p2 = static_cast<RooAbsReal *>(formula->getParameter(token2.c_str()));
778 RooAbsReal *p3 = static_cast<RooAbsReal *>(formula->getParameter(token3.c_str()));
779 RooRealVar *v2 = dynamic_cast<RooRealVar *>(p2);
780 RooRealVar *v3 = dynamic_cast<RooRealVar *>(p3);
781
782 auto *constr2 = findConstraint(v2);
783 auto *constr3 = findConstraint(v3);
784
785 if (constr2 && !p3) {
786 sys.name = p2->GetName();
787 sys.param = p2;
788 sys.high = sign * std::stod(token3);
789 sys.low = -sign * std::stod(token3);
790 } else if (!p2 && constr3) {
791 sys.name = p3->GetName();
792 sys.param = p3;
793 sys.high = sign * std::stod(token2);
794 sys.low = -sign * std::stod(token2);
795 } else if (constr2 && p3 && !constr3) {
796 sys.name = v2->GetName();
797 sys.param = v2;
798 sys.high = sign * p3->getVal();
799 sys.low = -sign * p3->getVal();
800 } else if (p2 && !constr2 && constr3) {
801 sys.name = v3->GetName();
802 sys.param = v3;
803 sys.high = sign * p2->getVal();
804 sys.low = -sign * p2->getVal();
805 }
806
807 // interpolation code 1 means linear, which is what we have here
808 sys.interpolationCode = 1;
809
810 erasePrefix(sys.name, "alpha_");
811 }
812 return sys;
813}
814
815void collectElements(RooArgSet &elems, RooAbsArg *arg)
816{
817 if (auto prod = dynamic_cast<RooProduct *>(arg)) {
818 for (const auto &e : prod->components()) {
819 collectElements(elems, e);
820 }
821 } else {
822 elems.add(*arg);
823 }
824}
825
826struct Sample {
827 std::string name;
828 std::vector<double> hist;
829 std::vector<double> histError;
830 std::vector<NormFactor> normfactors;
831 std::vector<NormSys> normsys;
832 std::vector<HistoSys> histosys;
833 std::vector<ShapeSys> shapesys;
834 std::vector<GenericElement> tmpElements;
835 std::vector<GenericElement> otherElements;
836 bool useBarlowBeestonLight = false;
837 std::vector<RooAbsReal *> staterrorParameters;
838 TClass *barlowBeestonLightConstraintType = RooPoisson::Class();
839 Sample(const std::string &n) : name{n} {}
840};
841
842void addNormFactor(RooRealVar const *par, Sample &sample, RooWorkspace *ws)
843{
844 std::string parname = par->GetName();
845 bool isConstrained = false;
846 for (RooAbsArg const *pdf : ws->allPdfs()) {
847 if (auto gauss = dynamic_cast<RooGaussian const *>(pdf)) {
848 if (parname == gauss->getX().GetName()) {
849 sample.normfactors.emplace_back(*par, gauss);
850 isConstrained = true;
851 }
852 }
853 }
854 if (!isConstrained)
855 sample.normfactors.emplace_back(*par);
856}
857
858namespace {
859
860bool verbose = false;
861
862}
863
864struct Channel {
865 std::string name;
866 std::vector<Sample> samples;
867 std::map<int, double> tot_yield;
868 std::map<int, double> tot_yield2;
869 std::map<int, double> rel_errors;
870 RooArgSet const *varSet = nullptr;
871 long unsigned int nBins = 0;
872};
873
875{
876 Channel channel;
877
878 RooWorkspace *ws = tool->workspace();
879
880 channel.name = pdfname;
881 erasePrefix(channel.name, "model_");
882 eraseSuffix(channel.name, "_model");
883
884 for (size_t sampleidx = 0; sampleidx < sumpdf->funcList().size(); ++sampleidx) {
885 PiecewiseInterpolation *pip = nullptr;
886 std::vector<ParamHistFunc *> phfs;
887
888 const auto func = sumpdf->funcList().at(sampleidx);
889 Sample sample(func->GetName());
890 erasePrefix(sample.name, "L_x_");
891 eraseSuffix(sample.name, "_shapes");
892 eraseSuffix(sample.name, "_" + channel.name);
893 erasePrefix(sample.name, pdfname + "_");
894
895 auto updateObservables = [&](RooDataHist const &dataHist) {
896 if (channel.varSet == nullptr) {
897 channel.varSet = dataHist.get();
898 channel.nBins = dataHist.numEntries();
899 }
900 if (sample.hist.empty()) {
901 auto *w = dataHist.weightArray();
902 sample.hist.assign(w, w + dataHist.numEntries());
903 }
904 };
905 auto processElements = [&](const auto &elements, auto &&self) -> void {
906 for (RooAbsArg *e : elements) {
907 if (TString(e->GetName()).Contains("binWidth")) {
908 // The bin width modifiers are handled separately. We can't just
909 // check for the RooBinWidthFunction type here, because prior to
910 // ROOT 6.26, the multiplication with the inverse bin width was
911 // done in a different way (like a normfactor with a RooRealVar,
912 // but it was stored in the dataset).
913 // Fortunately, the name was similar, so we can match the modifier
914 // name.
915 } else if (auto constVar = dynamic_cast<RooConstVar *>(e)) {
916 if (constVar->getVal() != 1.) {
917 sample.normfactors.emplace_back(*constVar);
918 }
919 } else if (auto par = dynamic_cast<RooRealVar *>(e)) {
920 addNormFactor(par, sample, ws);
921 } else if (auto hf = dynamic_cast<const RooHistFunc *>(e)) {
922 updateObservables(hf->dataHist());
923 } else if (auto phf = dynamic_cast<ParamHistFunc *>(e)) {
924 phfs.push_back(phf);
925 } else if (auto fip = dynamic_cast<RooStats::HistFactory::FlexibleInterpVar *>(e)) {
926 // some (modified) histfactory models have several instances of FlexibleInterpVar
927 // we collect and merge them
928 for (size_t i = 0; i < fip->variables().size(); ++i) {
929 RooAbsReal *var = static_cast<RooAbsReal *>(fip->variables().at(i));
930 std::string sysname(var->GetName());
931 erasePrefix(sysname, "alpha_");
932 const auto *constraint = findConstraint(var);
933 if (!constraint && !var->isConstant()) {
934 RooJSONFactoryWSTool::error("cannot find constraint for " + std::string(var->GetName()));
935 } else {
936 sample.normsys.emplace_back(sysname, var, fip->high()[i], fip->low()[i],
937 fip->interpolationCodes()[i], constraint);
938 }
939 }
940 } else if (!pip && (pip = dynamic_cast<PiecewiseInterpolation *>(e))) {
941 // nothing to do here, already assigned
942 } else if (RooFormulaVar *formula = dynamic_cast<RooFormulaVar *>(e)) {
943 // people do a lot of fancy stuff with RooFormulaVar, like including NormSys via explicit formulae.
944 // let's try to decompose it into building blocks
945 TString expression(formula->expression());
946 for (size_t i = formula->nParameters(); i--;) {
947 const RooAbsArg *p = formula->getParameter(i);
948 expression.ReplaceAll(("x[" + std::to_string(i) + "]").c_str(), p->GetName());
949 expression.ReplaceAll(("@" + std::to_string(i)).c_str(), p->GetName());
950 }
951 auto components = splitTopLevelProduct(expression.Data());
952 if (components.size() == 0) {
953 // it's not a product, let's just treat it as an unknown element
954 sample.otherElements.push_back(formula);
955 } else {
956 // it is a prododuct, we can try to handle the elements separately
957 std::vector<RooAbsArg *> realComponents;
958 int idx = 0;
959 for (auto &comp : components) {
960 // check if this is a trivial element of a product, we can treat it as its own modifier
961 auto *part = formula->getParameter(comp.c_str());
962 if (part) {
963 realComponents.push_back(part);
964 continue;
965 }
966 // check if this is an attempt at explicitly encoding an overallSys
967 auto normsys = parseOverallModifierFormula(comp, formula);
968 if (normsys.param) {
969 sample.normsys.emplace_back(std::move(normsys));
970 continue;
971 }
972
973 // this is something non-trivial, let's deal with it separately
974 std::string name = std::string(formula->GetName()) + "_part" + std::to_string(idx);
975 ++idx;
976 auto *var = new RooFormulaVar(name.c_str(), name.c_str(), comp.c_str(), formula->dependents());
977 sample.tmpElements.push_back({var});
978 }
979 self(realComponents, self);
980 }
981 } else if (auto real = dynamic_cast<RooAbsReal *>(e)) {
982 sample.otherElements.push_back(real);
983 }
984 }
985 };
986
987 RooArgSet elems;
988 collectElements(elems, func);
989 collectElements(elems, sumpdf->coefList().at(sampleidx));
991
992 // see if we can get the observables
993 if (pip) {
994 if (auto nh = dynamic_cast<RooHistFunc const *>(pip->nominalHist())) {
995 updateObservables(nh->dataHist());
996 }
997 }
998
999 // sort and configure norms
1000 sortByName(sample.normfactors);
1001 sortByName(sample.normsys);
1002
1003 // sort and configure the histosys
1004 if (pip) {
1005 for (size_t i = 0; i < pip->paramList().size(); ++i) {
1006 RooAbsReal *var = static_cast<RooAbsReal *>(pip->paramList().at(i));
1007 std::string sysname(var->GetName());
1008 erasePrefix(sysname, "alpha_");
1009 if (auto lo = dynamic_cast<RooHistFunc *>(pip->lowList().at(i))) {
1010 if (auto hi = dynamic_cast<RooHistFunc *>(pip->highList().at(i))) {
1011 const auto *constraint = findConstraint(var);
1012 if (!constraint && !var->isConstant()) {
1013 RooJSONFactoryWSTool::error("cannot find constraint for " + std::string(var->GetName()));
1014 } else {
1015 sample.histosys.emplace_back(sysname, var, lo, hi, constraint);
1016 }
1017 }
1018 }
1019 }
1020 sortByName(sample.histosys);
1021 }
1022
1023 for (ParamHistFunc *phf : phfs) {
1024 if (startsWith(std::string(phf->GetName()), "mc_stat_")) { // MC stat uncertainty
1025 int idx = 0;
1026 for (const auto &g : phf->paramList()) {
1027 sample.staterrorParameters.push_back(static_cast<RooRealVar *>(g));
1028 ++idx;
1029 RooAbsPdf *constraint = findConstraint(g);
1030 if (channel.tot_yield.find(idx) == channel.tot_yield.end()) {
1031 channel.tot_yield[idx] = 0;
1032 channel.tot_yield2[idx] = 0;
1033 }
1034 channel.tot_yield[idx] += sample.hist[idx - 1];
1035 channel.tot_yield2[idx] += (sample.hist[idx - 1] * sample.hist[idx - 1]);
1036 if (constraint) {
1037 sample.barlowBeestonLightConstraintType = constraint->IsA();
1038 if (RooPoisson *constraint_p = dynamic_cast<RooPoisson *>(constraint)) {
1039 double erel = 1. / std::sqrt(constraint_p->getX().getVal());
1040 channel.rel_errors[idx] = erel;
1041 } else if (RooGaussian *constraint_g = dynamic_cast<RooGaussian *>(constraint)) {
1042 double erel = constraint_g->getSigma().getVal() / constraint_g->getMean().getVal();
1043 channel.rel_errors[idx] = erel;
1044 } else {
1046 "currently, only RooPoisson and RooGaussian are supported as constraint types");
1047 }
1048 }
1049 }
1050 sample.useBarlowBeestonLight = true;
1051 } else { // other ShapeSys
1052 ShapeSys sys(phf->GetName());
1053 erasePrefix(sys.name, channel.name + "_");
1054 bool isshapesys = eraseSuffix(sys.name, "_ShapeSys") || eraseSuffix(sys.name, "_shapeSys");
1055 bool isshapefactor = eraseSuffix(sys.name, "_ShapeFactor") || eraseSuffix(sys.name, "_shapeFactor");
1056
1057 for (const auto &g : phf->paramList()) {
1058 sys.parameters.push_back(static_cast<RooRealVar *>(g));
1059 RooAbsPdf *constraint = nullptr;
1060 if (isshapesys) {
1061 constraint = findConstraint(g);
1062 if (!constraint)
1063 constraint = ws->pdf(constraintName(g->GetName()));
1064 if (!constraint && !g->isConstant()) {
1065 RooJSONFactoryWSTool::error("cannot find constraint for " + std::string(g->GetName()));
1066 }
1067 } else if (!isshapefactor) {
1068 RooJSONFactoryWSTool::error("unknown type of shapesys " + std::string(phf->GetName()));
1069 }
1070 if (!constraint) {
1071 sys.constraints.push_back(0.0);
1072 } else if (auto constraint_p = dynamic_cast<RooPoisson *>(constraint)) {
1073 sys.constraints.push_back(1. / std::sqrt(constraint_p->getX().getVal()));
1074 if (!sys.constraint) {
1075 sys.constraintType = RooPoisson::Class();
1076 }
1077 } else if (auto constraint_g = dynamic_cast<RooGaussian *>(constraint)) {
1078 sys.constraints.push_back(constraint_g->getSigma().getVal() / constraint_g->getMean().getVal());
1079 if (!sys.constraint) {
1080 sys.constraintType = RooGaussian::Class();
1081 }
1082 }
1083 }
1084 sample.shapesys.emplace_back(std::move(sys));
1085 }
1086 }
1087 sortByName(sample.shapesys);
1088
1089 // add the sample
1090 channel.samples.emplace_back(std::move(sample));
1091 }
1092
1093 sortByName(channel.samples);
1094 return channel;
1095}
1096
1097void configureStatError(Channel &channel)
1098{
1099 for (auto &sample : channel.samples) {
1100 if (sample.useBarlowBeestonLight) {
1101 sample.histError.resize(sample.hist.size());
1102 for (auto bin : channel.rel_errors) {
1103 // reverse engineering the correct partial error
1104 // the (arbitrary) convention used here is that all samples should have the same relative error
1105 const int i = bin.first;
1106 const double relerr_tot = bin.second;
1107 const double count = sample.hist[i - 1];
1108 // this reconstruction is inherently imprecise, so we truncate it at some decimal places to make sure that
1109 // we don't carry around too many useless digits
1110 sample.histError[i - 1] =
1111 round_prec(relerr_tot * channel.tot_yield[i] / std::sqrt(channel.tot_yield2[i]) * count, 7);
1112 }
1113 }
1114 }
1115}
1116
1118{
1119 bool observablesWritten = false;
1120 for (const auto &sample : channel.samples) {
1121
1122 elem["type"] << "histfactory_dist";
1123
1124 auto &s = RooJSONFactoryWSTool::appendNamedChild(elem["samples"], sample.name);
1125
1126 auto &modifiers = s["modifiers"];
1127 modifiers.set_seq();
1128
1129 for (const auto &nf : sample.normfactors) {
1130 auto &mod = modifiers.append_child();
1131 mod.set_map();
1132 mod["name"] << nf.name;
1133 mod["parameter"] << nf.param->GetName();
1134 mod["type"] << "normfactor";
1135 if (nf.constraint) {
1136 mod["constraint_name"] << nf.constraint->GetName();
1137 tool->queueExport(*nf.constraint);
1138 }
1139 }
1140
1141 for (const auto &sys : sample.normsys) {
1142 auto &mod = modifiers.append_child();
1143 mod.set_map();
1144 mod["name"] << sys.name;
1145 mod["type"] << "normsys";
1146 mod["parameter"] << sys.param->GetName();
1147 if (sys.interpolationCode != 4) {
1148 mod["interpolation"] << sys.interpolationCode;
1149 }
1150 if (sys.constraint) {
1151 mod["constraint_name"] << sys.constraint->GetName();
1152 } else if (sys.constraintType) {
1153 mod["constraint_type"] << toString(sys.constraintType);
1154 }
1155 auto &data = mod["data"].set_map();
1156 data["lo"] << sys.low;
1157 data["hi"] << sys.high;
1158 }
1159
1160 for (const auto &sys : sample.histosys) {
1161 auto &mod = modifiers.append_child();
1162 mod.set_map();
1163 mod["name"] << sys.name;
1164 mod["type"] << "histosys";
1165 mod["parameter"] << sys.param->GetName();
1166 if (sys.constraint) {
1167 mod["constraint_name"] << sys.constraint->GetName();
1168 } else if (sys.constraintType) {
1169 mod["constraint_type"] << toString(sys.constraintType);
1170 }
1171 auto &data = mod["data"].set_map();
1172 if (channel.nBins != sys.low.size() || channel.nBins != sys.high.size()) {
1173 std::stringstream ss;
1174 ss << "inconsistent binning: " << channel.nBins << " bins expected, but " << sys.low.size() << "/"
1175 << sys.high.size() << " found in nominal histogram errors!";
1176 RooJSONFactoryWSTool::error(ss.str().c_str());
1177 }
1178 RooJSONFactoryWSTool::exportArray(channel.nBins, sys.low.data(), data["lo"].set_map()["contents"]);
1179 RooJSONFactoryWSTool::exportArray(channel.nBins, sys.high.data(), data["hi"].set_map()["contents"]);
1180 }
1181
1182 for (const auto &sys : sample.shapesys) {
1183 auto &mod = modifiers.append_child();
1184 mod.set_map();
1185 mod["name"] << sys.name;
1186 mod["type"] << "shapesys";
1187 optionallyExportGammaParameters(mod, sys.name, sys.parameters);
1188 if (sys.constraint) {
1189 mod["constraint_name"] << sys.constraint->GetName();
1190 } else if (sys.constraintType) {
1191 mod["constraint_type"] << toString(sys.constraintType);
1192 }
1193 if (sys.constraint || sys.constraintType) {
1194 auto &vals = mod["data"].set_map()["vals"];
1195 vals.fill_seq(sys.constraints);
1196 } else {
1197 auto &vals = mod["data"].set_map()["vals"];
1198 vals.set_seq();
1199 for (std::size_t i = 0; i < sys.parameters.size(); ++i) {
1200 vals.append_child() << 0;
1201 }
1202 }
1203 }
1204
1205 for (const auto &other : sample.otherElements) {
1206 auto &mod = modifiers.append_child();
1207 mod.set_map();
1208 mod["name"] << other.name;
1209 mod["type"] << "custom";
1210 }
1211 for (const auto &other : sample.tmpElements) {
1212 auto &mod = modifiers.append_child();
1213 mod.set_map();
1214 mod["name"] << other.name;
1215 mod["type"] << "custom";
1216 }
1217
1218 if (sample.useBarlowBeestonLight) {
1219 auto &mod = modifiers.append_child();
1220 mod.set_map();
1221 mod["name"] << ::Literals::staterror;
1222 mod["type"] << ::Literals::staterror;
1223 optionallyExportGammaParameters(mod, "stat_" + channel.name, sample.staterrorParameters);
1224 mod["constraint_type"] << toString(sample.barlowBeestonLightConstraintType);
1225 }
1226
1227 if (!observablesWritten) {
1228 auto &output = elem["axes"].set_seq();
1229 for (auto *obs : static_range_cast<RooRealVar *>(*channel.varSet)) {
1230 auto &output.append_child().set_map();
1231 std::string name = obs->GetName();
1233 name;
1234 writeAxis(out, *obs);
1235 }
1236 observablesWritten = true;
1237 }
1238 auto &dataNode = s["data"].set_map();
1239 if (channel.nBins != sample.hist.size()) {
1240 std::stringstream ss;
1241 ss << "inconsistent binning: " << channel.nBins << " bins expected, but " << sample.hist.size()
1242 << " found in nominal histogram!";
1243 RooJSONFactoryWSTool::error(ss.str().c_str());
1244 }
1245 RooJSONFactoryWSTool::exportArray(channel.nBins, sample.hist.data(), dataNode["contents"]);
1246 if (!sample.histError.empty()) {
1247 if (channel.nBins != sample.histError.size()) {
1248 std::stringstream ss;
1249 ss << "inconsistent binning: " << channel.nBins << " bins expected, but " << sample.histError.size()
1250 << " found in nominal histogram errors!";
1251 RooJSONFactoryWSTool::error(ss.str().c_str());
1252 }
1253 RooJSONFactoryWSTool::exportArray(channel.nBins, sample.histError.data(), dataNode["errors"]);
1254 }
1255 }
1256
1257 return true;
1258}
1259
1260std::vector<RooAbsPdf *> findLostConstraints(const Channel &channel, const std::vector<RooAbsPdf *> &constraints)
1261{
1262 // collect all the vars that are used by the model
1263 std::set<const RooAbsReal *> vars;
1264 for (const auto &sample : channel.samples) {
1265 for (const auto &nf : sample.normfactors) {
1266 vars.insert(nf.param);
1267 }
1268 for (const auto &sys : sample.normsys) {
1269 vars.insert(sys.param);
1270 }
1271
1272 for (const auto &sys : sample.histosys) {
1273 vars.insert(sys.param);
1274 }
1275 for (const auto &sys : sample.shapesys) {
1276 for (const auto &par : sys.parameters) {
1277 vars.insert(par);
1278 }
1279 }
1280 if (sample.useBarlowBeestonLight) {
1281 for (const auto &par : sample.staterrorParameters) {
1282 vars.insert(par);
1283 }
1284 }
1285 }
1286
1287 // check if there is any constraint present that is unrelated to these vars
1288 std::vector<RooAbsPdf *> lostConstraints;
1289 for (auto *pdf : constraints) {
1290 bool related = false;
1291 for (const auto *var : vars) {
1292 if (pdf->dependsOn(*var)) {
1293 related = true;
1294 }
1295 }
1296 if (!related) {
1297 lostConstraints.push_back(pdf);
1298 }
1299 }
1300 // return the constraints that would be "lost" when exporting the model
1301 return lostConstraints;
1302}
1303
1305 std::vector<RooAbsPdf *> constraints, JSONNode &elem)
1306{
1307 // some preliminary checks
1308 if (!sumpdf) {
1309 if (verbose) {
1310 std::cout << pdfname << " is not a sumpdf" << std::endl;
1311 }
1312 return false;
1313 }
1314
1315 for (RooAbsArg *sample : sumpdf->funcList()) {
1316 if (!dynamic_cast<RooProduct *>(sample) && !dynamic_cast<RooRealSumPdf *>(sample)) {
1317 if (verbose)
1318 std::cout << "sample " << sample->GetName() << " is no RooProduct or RooRealSumPdf in " << pdfname
1319 << std::endl;
1320 return false;
1321 }
1322 }
1323
1324 auto channel = readChannel(tool, pdfname, sumpdf);
1325
1326 // sanity checks
1327 if (channel.samples.size() == 0)
1328 return false;
1329 for (auto &sample : channel.samples) {
1330 if (sample.hist.empty()) {
1331 return false;
1332 }
1333 }
1334
1335 // stat error handling
1336 configureStatError(channel);
1337
1338 auto lostConstraints = findLostConstraints(channel, constraints);
1339 // Export all the lost constraints
1340 for (const auto *constraint : lostConstraints) {
1342 "losing constraint term '" + std::string(constraint->GetName()) +
1343 "', implicit constraints are not supported by HS3 yet! The term will appear in the HS3 file, but will not be "
1344 "picked up when creating a likelihood from it! You will have to add it manually as an external constraint.");
1345 tool->queueExport(*constraint);
1346 }
1347
1348 // Export all the regular modifiers
1349 for (const auto &sample : channel.samples) {
1350 for (auto &modifier : sample.normfactors) {
1351 if (modifier.constraint) {
1352 tool->queueExport(*modifier.constraint);
1353 }
1354 }
1355 for (auto &modifier : sample.normsys) {
1356 if (modifier.constraint) {
1357 tool->queueExport(*modifier.constraint);
1358 }
1359 }
1360 for (auto &modifier : sample.histosys) {
1361 if (modifier.constraint) {
1362 tool->queueExport(*modifier.constraint);
1363 }
1364 }
1365 }
1366
1367 // Export all the custom modifiers
1368 for (const auto &sample : channel.samples) {
1369 for (auto &modifier : sample.otherElements) {
1370 tool->queueExport(*modifier.function);
1371 }
1372 for (auto &modifier : sample.tmpElements) {
1373 tool->queueExportTemporary(modifier.function);
1374 }
1375 }
1376
1377 // Export all model parameters
1378 RooArgSet parameters;
1379 sumpdf->getParameters(channel.varSet, parameters);
1380 for (RooAbsArg *param : parameters) {
1381 // This should exclude the global observables
1382 if (!startsWith(std::string{param->GetName()}, "nom_")) {
1383 tool->queueExport(*param);
1384 }
1385 }
1386
1387 return exportChannel(tool, channel, elem);
1388}
1389
1390class HistFactoryStreamer_ProdPdf : public RooFit::JSONIO::Exporter {
1391public:
1392 bool autoExportDependants() const override { return false; }
1394 {
1395 std::vector<RooAbsPdf *> constraints;
1396 RooRealSumPdf *sumpdf = nullptr;
1397 for (RooAbsArg *v : prodpdf->pdfList()) {
1398 RooAbsPdf *pdf = static_cast<RooAbsPdf *>(v);
1399 auto thispdf = dynamic_cast<RooRealSumPdf *>(pdf);
1400 if (thispdf) {
1401 if (!sumpdf)
1402 sumpdf = thispdf;
1403 else
1404 return false;
1405 } else {
1406 constraints.push_back(pdf);
1407 }
1408 }
1409 if (!sumpdf)
1410 return false;
1411
1412 bool ok = tryExportHistFactory(tool, prodpdf->GetName(), sumpdf, constraints, elem);
1413 return ok;
1414 }
1415 std::string const &key() const override
1416 {
1417 static const std::string keystring = "histfactory_dist";
1418 return keystring;
1419 }
1420 bool exportObject(RooJSONFactoryWSTool *tool, const RooAbsArg *p, JSONNode &elem) const override
1421 {
1422 return tryExport(tool, static_cast<const RooProdPdf *>(p), elem);
1423 }
1424};
1425
1426class HistFactoryStreamer_SumPdf : public RooFit::JSONIO::Exporter {
1427public:
1428 bool autoExportDependants() const override { return false; }
1430 {
1431 std::vector<RooAbsPdf *> constraints;
1432 return tryExportHistFactory(tool, sumpdf->GetName(), sumpdf, constraints, elem);
1433 }
1434 std::string const &key() const override
1435 {
1436 static const std::string keystring = "histfactory_dist";
1437 return keystring;
1438 }
1439 bool exportObject(RooJSONFactoryWSTool *tool, const RooAbsArg *p, JSONNode &elem) const override
1440 {
1441 return tryExport(tool, static_cast<const RooRealSumPdf *>(p), elem);
1442 }
1443};
1444
1445STATIC_EXECUTE([]() {
1446 using namespace RooFit::JSONIO;
1447
1448 registerImporter<HistFactoryImporter>("histfactory_dist", true);
1450 registerImporter<FlexibleInterpVarFactory>("interpolation0d", true);
1455});
1456
1457} // namespace
bool startsWith(std::string_view str, std::string_view prefix)
bool endsWith(std::string_view str, std::string_view suffix)
#define d(i)
Definition RSha256.hxx:102
#define c(i)
Definition RSha256.hxx:101
#define g(i)
Definition RSha256.hxx:105
#define h(i)
Definition RSha256.hxx:106
#define e(i)
Definition RSha256.hxx:103
ROOT::RRangeCast< T, false, Range_t > static_range_cast(Range_t &&coll)
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void funcs
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t modifier
char name[80]
Definition TGX11.cxx:110
#define hi
TClass * IsA() const override
Definition TStringLong.h:20
A class which maps the current values of a RooRealVar (or a set of RooRealVars) to one of a number of...
The PiecewiseInterpolation is a class that can morph distributions into each other,...
static TClass * Class()
void setPositiveDefinite(bool flag=true)
const_iterator begin() const
const_iterator end() const
Common abstract base class for objects that represent a value and a "shape" in RooFit.
Definition RooAbsArg.h:77
bool dependsOn(const RooAbsCollection &serverList, const RooAbsArg *ignoreArg=nullptr, bool valueOnly=false) const
Test whether we depend on (ie, are served by) any object in the specified collection.
bool isConstant() const
Check if the "Constant" attribute is set.
Definition RooAbsArg.h:304
Abstract container object that can hold multiple RooAbsArg objects.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
TClass * IsA() const override
Definition RooAbsPdf.h:351
Int_t numBins(const char *rangeName=nullptr) const override
virtual double getMax(const char *name=nullptr) const
Get maximum of currently defined range.
virtual double getMin(const char *name=nullptr) const
Get minimum of currently defined range.
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:24
Returns the bin width (or volume) given a RooHistFunc.
Represents a constant real-valued object.
Definition RooConstVar.h:23
Container class to hold N-dimensional binned data.
Definition RooDataHist.h:40
virtual JSONNode & set_seq()=0
A RooFormulaVar is a generic implementation of a real-valued object, which takes a RooArgList of serv...
RooAbsArg * getParameter(const char *name) const
Return pointer to parameter with given name.
const char * expression() const
const RooArgList & dependents() const
size_t nParameters() const
Return the number of parameters.
Plain Gaussian p.d.f.
Definition RooGaussian.h:24
static TClass * Class()
A real-valued function sampled from a multidimensional histogram.
Definition RooHistFunc.h:31
When using RooFit, statistical models can be conveniently handled and stored as a RooWorkspace.
static void fillSeq(RooFit::Detail::JSONNode &node, RooAbsCollection const &coll, size_t nMax=-1)
static std::unique_ptr< RooDataHist > readBinnedData(const RooFit::Detail::JSONNode &n, const std::string &namecomp, RooArgSet const &vars)
Read binned data from the JSONNode and create a RooDataHist object.
static RooFit::Detail::JSONNode & appendNamedChild(RooFit::Detail::JSONNode &node, std::string const &name)
static void exportArray(std::size_t n, double const *contents, RooFit::Detail::JSONNode &output)
Export an array of doubles to a JSONNode.
static bool testValidName(const std::string &str, bool forcError)
static void error(const char *s)
Writes an error message to the RooFit message service and throws a runtime_error.
static std::string name(const RooFit::Detail::JSONNode &n)
static std::ostream & warning(const std::string &s)
Writes a warning message to the RooFit message service.
static RooArgSet readAxes(const RooFit::Detail::JSONNode &node)
Read axes from the JSONNode and create a RooArgSet representing them.
RooFit Lognormal PDF.
static TClass * Class()
Poisson pdf.
Definition RooPoisson.h:19
static TClass * Class()
Efficient implementation of a product of PDFs of the form.
Definition RooProdPdf.h:36
static TClass * Class()
Represents the product of a given set of RooAbsReal objects.
Definition RooProduct.h:29
Implements a PDF constructed from a sum of functions:
static TClass * Class()
Variable that can be changed from the outside.
Definition RooRealVar.h:37
void setError(double value)
Definition RooRealVar.h:60
const RooAbsBinning & getBinning(const char *name=nullptr, bool verbose=true, bool createOnTheFly=false) const override
Return binning definition with name.
This class encapsulates all information for the statistical interpretation of one experiment.
Definition Channel.h:30
Configuration for a constrained, coherent shape variation of affected samples.
Configuration for an un- constrained overall systematic to scale sample normalisations.
Definition Systematics.h:63
Constrained bin-by-bin variation of affected histogram.
Persistable container for RooFit projects.
RooAbsPdf * pdf(RooStringView name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
RooAbsReal * function(RooStringView name) const
Retrieve function (RooAbsReal) with given name. Note that all RooAbsPdfs are also RooAbsReals....
RooFactoryWSTool & factory()
Return instance to factory tool.
bool import(const RooAbsArg &arg, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}, const RooCmdArg &arg9={})
Import a RooAbsArg object, e.g.
TClass instances represent classes, structs and namespaces in the ROOT type system.
Definition TClass.h:84
const char * GetName() const override
Returns name of object.
Definition TNamed.h:49
Basic string class.
Definition TString.h:139
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:632
RooCmdArg RecycleConflictNodes(bool flag=true)
RooCmdArg Conditional(const RooArgSet &pdfSet, const RooArgSet &depSet, bool depsAreCond=false)
const Int_t n
Definition legend1.C:16
double gamma(double x)
void function(const Char_t *name_, T fun, const Char_t *docstring=0)
Definition RExports.h:167
CreateGammaConstraintsOutput createGammaConstraints(RooArgList const &paramList, std::span< const double > relSigmas, double minSigma, Constraint::Type type)
#define STATIC_EXECUTE(MY_FUNC)
TLine l
Definition textangle.C:4
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345
static void output()