Skip to content

Commit ec5952f

Browse files
committed
[HS3] Fix HistFactory MC stat uncertainty export
1 parent 9916067 commit ec5952f

1 file changed

Lines changed: 33 additions & 2 deletions

File tree

roofit/hs3/src/JSONFactories_HistFactory.cxx

Lines changed: 33 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -283,6 +283,36 @@ getOrCreateConstraint(RooJSONFactoryWSTool &tool, const JSONNode &mod, RooRealVa
283283
"'");
284284
}
285285
}
286+
double poissonTau(RooPoisson const &constraint, RooAbsArg const &gamma)
287+
{
288+
auto const *mean = dynamic_cast<RooProduct const *>(&constraint.getMean());
289+
if (!mean) {
290+
RooJSONFactoryWSTool::error("Poisson gamma constraint mean is not a RooProduct: " +
291+
std::string(constraint.GetName()));
292+
}
293+
294+
for (RooAbsArg *arg : mean->servers()) {
295+
if (arg == &gamma) {
296+
continue;
297+
}
298+
299+
if (auto const *tau = dynamic_cast<RooConstVar const *>(arg)) {
300+
return tau->getVal();
301+
}
302+
303+
// Imported workspaces can sometimes represent
304+
// constants as constant RooRealVars.
305+
if (auto const *real = dynamic_cast<RooAbsReal const *>(arg)) {
306+
if (real->isConstant() || endsWith(std::string(real->GetName()), "_tau")) {
307+
return real->getVal();
308+
}
309+
}
310+
}
311+
312+
RooJSONFactoryWSTool::error("Could not find tau component in Poisson gamma constraint mean: " +
313+
std::string(constraint.GetName()));
314+
return std::numeric_limits<double>::quiet_NaN();
315+
}
286316

287317
bool importHistSample(RooJSONFactoryWSTool &tool, RooDataHist &dh, RooArgSet const &varlist,
288318
RooAbsArg const *mcStatObject, const std::string &fprefix, const JSONNode &p,
@@ -334,6 +364,7 @@ bool importHistSample(RooJSONFactoryWSTool &tool, RooDataHist &dh, RooArgSet con
334364
// this is dealt with at a different place, ignore it for now
335365
} else if (modtype == "normfactor") {
336366
RooRealVar &constrParam = getOrCreate<RooRealVar>(ws, sysname, 1., -3, 5);
367+
constrParam.setError(0.0);
337368
normElems.add(constrParam);
338369
if (mod.has_child("constraint_name") || mod.has_child("constraint_type")) {
339370
// for norm factors, constraints are optional
@@ -1060,7 +1091,7 @@ Channel readChannel(RooJSONFactoryWSTool *tool, const std::string &pdfname, cons
10601091
if (constraint) {
10611092
sample.barlowBeestonLightConstraintType = constraint->IsA();
10621093
if (RooPoisson *constraint_p = dynamic_cast<RooPoisson *>(constraint)) {
1063-
double erel = 1. / std::sqrt(constraint_p->getX().getVal());
1094+
double erel = 1. / std::sqrt(poissonTau(*constraint_p, *g));
10641095
channel.rel_errors[idx] = erel;
10651096
} else if (RooGaussian *constraint_g = dynamic_cast<RooGaussian *>(constraint)) {
10661097
double erel = constraint_g->getSigma().getVal() / constraint_g->getMean().getVal();
@@ -1094,7 +1125,7 @@ Channel readChannel(RooJSONFactoryWSTool *tool, const std::string &pdfname, cons
10941125
if (!constraint) {
10951126
sys.constraints.push_back(0.0);
10961127
} else if (auto constraint_p = dynamic_cast<RooPoisson *>(constraint)) {
1097-
sys.constraints.push_back(1. / std::sqrt(constraint_p->getX().getVal()));
1128+
sys.constraints.push_back(1. / std::sqrt(poissonTau(*constraint_p, *g)));
10981129
if (!sys.constraint) {
10991130
sys.constraintType = RooPoisson::Class();
11001131
}

0 commit comments

Comments
 (0)