51 if (binning.isUniform()) {
53 axis[
"min"] << obs.
getMin();
54 axis[
"max"] << obs.
getMax();
56 auto &edges = axis[
"edges"];
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;
71 int ndigits = std::floor(std::log10(std::abs(
d))) + 1 -
nSig;
73 if (std::abs(
d /
sf) < 2)
75 return sf * std::round(
d /
sf);
84void erasePrefix(std::string &str, std::string_view prefix)
87 str.erase(0, prefix.size());
94 str.erase(str.size() -
suffix.size());
104 std::sort(
coll.
begin(),
coll.
end(), [](
auto &
l,
auto &
r) { return l.name < r.name; });
110 for (
const auto &client :
gamma->clients()) {
111 if (
auto casted =
dynamic_cast<T *
>(client)) {
157 return "gamma_" +
sysname +
"_bin_" + std::to_string(i);
167 for (std::size_t i = 0; i < params.size(); ++i) {
168 std::string
name(params[i]->GetName());
182 nom.setConstant(
true);
194 const std::vector<std::string> &
parnames,
const std::vector<double> &vals,
201 for (std::size_t i = 0; i < vals.size(); ++i) {
208 if (constraintType !=
"Const") {
217 gamma->setConstant(
true);
226 if (!
comp.has_child(
"modifiers"))
228 for (
const auto &
mod :
comp[
"modifiers"].children()) {
229 if (
mod[
"type"].val() == ::Literals::staterror)
237 if (
comp.has_child(
"modifiers")) {
238 for (
const auto &
mod :
comp[
"modifiers"].children()) {
239 if (
mod[
"type"].val() == ::Literals::staterror)
244 ::Literals::staterror +
" modifier!");
261 param.
setError(gauss->getSigma().getVal());
265 std::cout <<
"creating new constraint for " << param << std::endl;
273 *
tool.workspace()->var(std::string(
"nom_") + param.
GetName()), 1.);
292 if (!
p.has_child(
"data")) {
308 if (
p.has_child(
"modifiers")) {
319 for (
const auto &
mod :
p[
"modifiers"].children()) {
320 std::string
const &
modtype =
mod[
"type"].val();
322 mod.has_child(
"name")
324 : (
mod.has_child(
"parameter") ?
mod[
"parameter"].val() :
"syst_" + std::to_string(idx));
328 }
else if (
modtype ==
"normfactor") {
331 if (
mod.has_child(
"constraint_name") ||
mod.has_child(
"constraint_type")) {
335 }
else if (
modtype ==
"normsys") {
343 if (
mod.has_child(
"interpolation")) {
346 double low =
data[
"lo"].val_double();
347 double high =
data[
"hi"].val_double();
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();
363 }
else if (
modtype ==
"histosys") {
377 }
else if (
modtype ==
"shapesys") {
380 std::vector<double> vals;
381 for (
const auto &
v :
mod[
"data"][
"vals"].children()) {
382 vals.push_back(
v.val_double());
385 for (
const auto &
v :
mod[
"parameters"].children()) {
391 std::string constraint(
mod.has_child(
"constraint_type") ?
mod[
"constraint_type"].val()
392 :
mod.has_child(
"constraint") ?
mod[
"constraint"].val()
396 }
else if (
modtype ==
"custom") {
420 v.setAllInterpCodes(4);
442 if (!
p.has_child(
"samples")) {
447 if (
p.has_child(::Literals::staterror)) {
448 auto &
staterr =
p[::Literals::staterror];
449 if (
staterr.has_child(
"relThreshold"))
451 if (
staterr.has_child(
"constraint_type"))
454 std::vector<double>
sumW;
455 std::vector<double>
sumW2;
461 std::vector<std::unique_ptr<RooDataHist>>
data;
462 for (
const auto &
comp :
p[
"samples"].children()) {
465 size_t nbins =
dh->numEntries();
472 for (
size_t i = 0; i < nbins; ++i) {
474 sumW2[i] +=
dh->weightSquared(i);
484 data.emplace_back(std::move(
dh));
493 std::vector<double>
errs(
sumW.size());
494 for (
size_t i = 0; i <
sumW.size(); ++i) {
508 for (
const auto &
comp :
p[
"samples"].children()) {
517 if (constraints.
empty()) {
532 std::string
const &key()
const override
534 static const std::string
keystring =
"interpolation0d";
540 elem[
"type"] << key();
541 elem[
"interpolationCodes"].fill_seq(
fip->interpolationCodes());
544 elem[
"high"].fill_seq(
fip->high(),
fip->variables().size());
545 elem[
"low"].fill_seq(
fip->low(),
fip->variables().size());
552 std::string
const &key()
const override
554 static const std::string
keystring =
"interpolation";
560 elem[
"type"] << key();
561 elem[
"interpolationCodes"].fill_seq(
pip->interpolationCodes());
562 elem[
"positiveDefinite"] <<
pip->positiveDefinite();
564 elem[
"nom"] <<
pip->nominalHist()->GetName();
583 pip.setPositiveDefinite(
p[
"positiveDefinite"].val_bool());
585 if (
p.has_child(
"interpolationCodes")) {
587 for (
auto const &node :
p[
"interpolationCodes"].children()) {
588 pip.setInterpCode(*
static_cast<RooAbsReal *
>(vars.at(i)), node.val_int(),
true);
602 if (!
p.has_child(
"high")) {
605 if (!
p.has_child(
"low")) {
608 if (!
p.has_child(
"nom")) {
612 double nom(
p[
"nom"].val_double());
616 std::vector<double> high;
619 std::vector<double> low;
622 if (vars.size() != low.size() || vars.size() != high.size()) {
624 "' has non-matching lengths of 'vars', 'high' and 'low'!");
629 if (
p.has_child(
"interpolationCodes")) {
631 for (
auto const &node :
p[
"interpolationCodes"].children()) {
632 fip.setInterpCode(*
static_cast<RooAbsReal *
>(vars.at(i)), node.val_int());
653 std::string
name =
"";
657 int interpolationCode = 4;
662 :
name(
n), param(
p), low(
l), high(
h), interpolationCode(i), constraint(
c), constraintType(
c->
IsA())
670 std::vector<double> low;
671 std::vector<double> high;
675 :
name(
n), param(
p), constraint(
c), constraintType(
c->
IsA())
677 low.assign(
l->dataHist().weightArray(),
l->dataHist().weightArray() +
l->dataHist().numEntries());
678 high.assign(
h->dataHist().weightArray(),
h->dataHist().weightArray() +
h->dataHist().numEntries());
683 std::vector<double> constraints;
684 std::vector<RooAbsReal *> parameters;
690struct GenericElement {
699 size_t end = s.size();
701 while (start < end && s[start] ==
'(' && s[end - 1] ==
')') {
704 for (
size_t i = start; i < end - 1; ++i) {
707 else if (s[i] ==
')')
709 if (
depth == 0 && i < end - 1) {
721 return s.substr(start, end - start);
726 std::vector<std::string>
parts;
731 for (
size_t i = 0; i <
expr.size(); ++i) {
735 }
else if (
c ==
')') {
737 }
else if (
c ==
'*' &&
depth == 0) {
739 std::string sub =
expr.substr(start, i - start);
749 std::string sub =
expr.substr(start);
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*$)");
769 if (std::regex_match(s, match, pattern)) {
770 if (match[1].str() ==
"-") {
774 std::string
token2 = match[2].str();
775 std::string
token3 = match[4].str();
786 sys.name =
p2->GetName();
791 sys.name =
p3->GetName();
796 sys.name =
v2->GetName();
798 sys.high =
sign *
p3->getVal();
799 sys.low = -
sign *
p3->getVal();
801 sys.name =
v3->GetName();
803 sys.high =
sign *
p2->getVal();
804 sys.low = -
sign *
p2->getVal();
808 sys.interpolationCode = 1;
817 if (
auto prod =
dynamic_cast<RooProduct *
>(arg)) {
818 for (
const auto &
e : prod->components()) {
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;
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);
855 sample.normfactors.emplace_back(*par);
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;
871 long unsigned int nBins = 0;
886 std::vector<ParamHistFunc *>
phfs;
896 if (channel.varSet ==
nullptr) {
897 channel.varSet = dataHist.get();
898 channel.nBins = dataHist.numEntries();
900 if (
sample.hist.empty()) {
901 auto *
w = dataHist.weightArray();
902 sample.hist.assign(
w,
w + dataHist.numEntries());
919 }
else if (
auto par =
dynamic_cast<RooRealVar *
>(
e)) {
928 for (
size_t i = 0; i <
fip->variables().
size(); ++i) {
937 fip->interpolationCodes()[i], constraint);
948 expression.ReplaceAll((
"x[" + std::to_string(i) +
"]").c_str(),
p->GetName());
949 expression.ReplaceAll((
"@" + std::to_string(i)).c_str(),
p->GetName());
952 if (components.size() == 0) {
954 sample.otherElements.push_back(formula);
957 std::vector<RooAbsArg *> realComponents;
959 for (
auto &
comp : components) {
963 realComponents.push_back(
part);
969 sample.normsys.emplace_back(std::move(normsys));
974 std::string
name = std::string(formula->
GetName()) +
"_part" + std::to_string(idx);
977 sample.tmpElements.push_back({var});
1005 for (
size_t i = 0; i <
pip->paramList().
size(); ++i) {
1009 if (
auto lo =
dynamic_cast<RooHistFunc *
>(
pip->lowList().at(i))) {
1026 for (
const auto &
g :
phf->paramList()) {
1030 if (channel.tot_yield.find(idx) == channel.tot_yield.end()) {
1031 channel.tot_yield[idx] = 0;
1032 channel.tot_yield2[idx] = 0;
1034 channel.tot_yield[idx] +=
sample.hist[idx - 1];
1035 channel.tot_yield2[idx] += (
sample.hist[idx - 1] *
sample.hist[idx - 1]);
1037 sample.barlowBeestonLightConstraintType = constraint->
IsA();
1040 channel.rel_errors[idx] =
erel;
1043 channel.rel_errors[idx] =
erel;
1046 "currently, only RooPoisson and RooGaussian are supported as constraint types");
1050 sample.useBarlowBeestonLight =
true;
1057 for (
const auto &
g :
phf->paramList()) {
1058 sys.parameters.push_back(
static_cast<RooRealVar *
>(
g));
1064 if (!constraint && !
g->isConstant()) {
1071 sys.constraints.push_back(0.0);
1073 sys.constraints.push_back(1. / std::sqrt(
constraint_p->getX().getVal()));
1074 if (!sys.constraint) {
1079 if (!sys.constraint) {
1084 sample.shapesys.emplace_back(std::move(sys));
1090 channel.samples.emplace_back(std::move(
sample));
1099 for (
auto &
sample : channel.samples) {
1100 if (
sample.useBarlowBeestonLight) {
1102 for (
auto bin : channel.rel_errors) {
1105 const int i = bin.first;
1107 const double count =
sample.hist[i - 1];
1110 sample.histError[i - 1] =
1120 for (
const auto &
sample : channel.samples) {
1122 elem[
"type"] <<
"histfactory_dist";
1129 for (
const auto &
nf :
sample.normfactors) {
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);
1141 for (
const auto &sys :
sample.normsys) {
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;
1150 if (sys.constraint) {
1151 mod[
"constraint_name"] << sys.constraint->GetName();
1152 }
else if (sys.constraintType) {
1153 mod[
"constraint_type"] << toString(sys.constraintType);
1155 auto &
data =
mod[
"data"].set_map();
1156 data[
"lo"] << sys.low;
1157 data[
"hi"] << sys.high;
1160 for (
const auto &sys :
sample.histosys) {
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);
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!";
1182 for (
const auto &sys :
sample.shapesys) {
1185 mod[
"name"] << sys.name;
1186 mod[
"type"] <<
"shapesys";
1188 if (sys.constraint) {
1189 mod[
"constraint_name"] << sys.constraint->GetName();
1190 }
else if (sys.constraintType) {
1191 mod[
"constraint_type"] << toString(sys.constraintType);
1193 if (sys.constraint || sys.constraintType) {
1194 auto &vals =
mod[
"data"].set_map()[
"vals"];
1195 vals.fill_seq(sys.constraints);
1197 auto &vals =
mod[
"data"].set_map()[
"vals"];
1199 for (std::size_t i = 0; i < sys.parameters.size(); ++i) {
1200 vals.append_child() << 0;
1209 mod[
"type"] <<
"custom";
1215 mod[
"type"] <<
"custom";
1218 if (
sample.useBarlowBeestonLight) {
1221 mod[
"name"] << ::Literals::staterror;
1222 mod[
"type"] << ::Literals::staterror;
1224 mod[
"constraint_type"] << toString(
sample.barlowBeestonLightConstraintType);
1230 auto &output.append_child().set_map();
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!";
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!";
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);
1268 for (
const auto &sys :
sample.normsys) {
1269 vars.insert(sys.param);
1272 for (
const auto &sys :
sample.histosys) {
1273 vars.insert(sys.param);
1275 for (
const auto &sys :
sample.shapesys) {
1276 for (
const auto &par : sys.parameters) {
1280 if (
sample.useBarlowBeestonLight) {
1281 for (
const auto &par :
sample.staterrorParameters) {
1289 for (
auto *pdf : constraints) {
1291 for (
const auto *var : vars) {
1292 if (pdf->dependsOn(*var)) {
1310 std::cout <<
pdfname <<
" is not a sumpdf" << std::endl;
1318 std::cout <<
"sample " <<
sample->GetName() <<
" is no RooProduct or RooRealSumPdf in " <<
pdfname
1327 if (channel.samples.size() == 0)
1329 for (
auto &
sample : channel.samples) {
1330 if (
sample.hist.empty()) {
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);
1349 for (
const auto &
sample : channel.samples) {
1368 for (
const auto &
sample : channel.samples) {
1379 sumpdf->getParameters(channel.varSet, parameters);
1383 tool->queueExport(*param);
1392 bool autoExportDependants()
const override {
return false; }
1395 std::vector<RooAbsPdf *> constraints;
1406 constraints.push_back(pdf);
1415 std::string
const &key()
const override
1417 static const std::string
keystring =
"histfactory_dist";
1428 bool autoExportDependants()
const override {
return false; }
1431 std::vector<RooAbsPdf *> constraints;
1434 std::string
const &key()
const override
1436 static const std::string
keystring =
"histfactory_dist";
bool startsWith(std::string_view str, std::string_view prefix)
bool endsWith(std::string_view str, std::string_view suffix)
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
TClass * IsA() const override
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,...
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.
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.
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.
TClass * IsA() const override
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...
RooArgList is a container object that can hold multiple RooAbsArg objects.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Returns the bin width (or volume) given a RooHistFunc.
Represents a constant real-valued object.
Container class to hold N-dimensional binned data.
virtual JSONNode & set_seq()=0
A real-valued function sampled from a multidimensional histogram.
Efficient implementation of a product of PDFs of the form.
Represents the product of a given set of RooAbsReal objects.
Implements a PDF constructed from a sum of functions:
Variable that can be changed from the outside.
void setError(double value)
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.
Configuration for a constrained, coherent shape variation of affected samples.
Configuration for an un- constrained overall systematic to scale sample normalisations.
std::string GetName() const
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.
const char * GetName() const override
Returns name of object.
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
RooCmdArg RecycleConflictNodes(bool flag=true)
RooCmdArg Conditional(const RooArgSet &pdfSet, const RooArgSet &depSet, bool depsAreCond=false)
void function(const Char_t *name_, T fun, const Char_t *docstring=0)
constexpr double defaultShapeSysGammaMax
constexpr double minShapeUncertainty
constexpr double defaultStatErrorGammaMax
constexpr double defaultGammaMin
CreateGammaConstraintsOutput createGammaConstraints(RooArgList const ¶mList, std::span< const double > relSigmas, double minSigma, Constraint::Type type)
#define STATIC_EXECUTE(MY_FUNC)
static uint64_t sum(uint64_t i)