Skip to content

Commit

Permalink
[RF][HF] Consistent limit and error configuration for gamma parameters
Browse files Browse the repository at this point in the history 
This commit unifyies the gamma parameter configuration in the original
HistFactory implementation and the new JSON serialization.
  • Loading branch information
@guitargeek
guitargeek committed Aug 31, 2023
1 parent 44efaba commit 3c68044
Show file tree
Hide file tree
Showing 5 changed files with 107 additions and 69 deletions.
14 changes: 14 additions & 0 deletions roofit/histfactory/inc/RooStats/HistFactory/Detail/HistFactoryImpl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,10 +16,22 @@
#include <RooGlobalFunc.h>
#include <RooWorkspace.h>

#include <ROOT/RSpan.hxx>

namespace RooStats {
namespace HistFactory {
namespace Detail {

namespace MagicConstants {

constexpr double defaultGammaMin = 0;
constexpr double defaultShapeFactorGammaMax = 1000;
constexpr double defaultShapeSysGammaMax = 10;
constexpr double defaultStatErrorGammaMax = 10;
constexpr double minShapeUncertainty = 0.0;

} // namespace MagicConstants

template <class Arg_t, class... Params_t>
Arg_t &getOrCreate(RooWorkspace &ws, std::string const &name, Params_t &&...params)
{
Expand All @@ -31,6 +43,8 @@ Arg_t &getOrCreate(RooWorkspace &ws, std::string const &name, Params_t &&...para
return *static_cast<Arg_t *>(ws.obj(name));
}

void configureConstrainedGammas(RooArgList const &gammas, std::span<double> relSigmas, double minSigma);

} // namespace Detail
} // namespace HistFactory
} // namespace RooStats
Expand Down
2 changes: 1 addition & 1 deletion roofit/histfactory/inc/RooStats/HistFactory/HistoToWorkspaceFactoryFast.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,7 @@ namespace RooStats{

TH1* MakeAbsolUncertaintyHist( const std::string& Name, const TH1* Hist );

RooArgList createStatConstraintTerms( RooWorkspace& proto,
RooArgList createGammaConstraintTerms( RooWorkspace& proto,
std::vector<std::string>& constraintTerms,
ParamHistFunc& paramHist, const TH1* uncertHist,
Constraint::Type type, double minSigma );
Expand Down
66 changes: 66 additions & 0 deletions roofit/histfactory/src/HistFactoryImpl.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,3 +11,69 @@
*/

#include <RooStats/HistFactory/Detail/HistFactoryImpl.h>

#include <RooMsgService.h>
#include <RooRealVar.h>

namespace RooStats {
namespace HistFactory {
namespace Detail {

/**
* \brief Configure constrained gamma parameters for fitting.
*
* This function configures constrained gamma parameters for fitting. If a
* given relative sigma is less than or equal to zero or below a threshold, the
* gamma parameter is set to be constant. The function also sets reasonable
* ranges for the gamma parameter and provides a reasonable starting point for
* pre-fit errors.
*
* @param gammas The gamma parameters to be configured.
* @param sigmaRel The relative sigma values to be used for configuring the
* limits and errors.
* @param minSigma The minimum relative sigma threshold. If a relative sigma is
* below this threshold, the gamma parameter is set to be
* constant.
*/
void configureConstrainedGammas(RooArgList const &gammas, std::span<double> relSigmas, double minSigma)
{
assert(gammas.size() == relSigmas.size());

for (std::size_t i = 0; i < gammas.size(); ++i) {
auto &gamma = *static_cast<RooRealVar *>(gammas.at(i));
double sigmaRel = relSigmas[i];

// If the sigma is zero, the parameter might as well be constant
if (sigmaRel <= 0) {
gamma.setConstant(true);
continue;
}

// Set reasonable ranges
gamma.setMax(1. + 5. * sigmaRel);
gamma.setMin(0.);

// Give reasonable starting point for pre-fit errors by setting it to the
// absolute sigma Mostly useful for pre-fit plotting.
// Note: in commit 2129c4d920 "[HF] Reduce verbosity of HistFactory."
// from 2020, there was a check added to do this only for Gaussian
// constrained parameters and for Poisson constrained parameters if they
// are stat errors without any justification. In the ROOT 6.30
// development cycle, this check got removed again to cause less surprise
// to the user.
gamma.setError(sigmaRel);

// If the sigma value is less than a supplied threshold, set the variable to
// constant
if (sigmaRel < minSigma) {
oocxcoutW(static_cast<TObject *>(nullptr), HistFactory)
<< "Warning: relative sigma " << sigmaRel << " for \"" << gamma.GetName() << "\" falls below threshold of "
<< minSigma << ". Setting: " << gamma.GetName() << " to constant" << std::endl;
gamma.setConstant(true);
}
}
}

} // namespace Detail
} // namespace HistFactory
} // namespace RooStats
67 changes: 18 additions & 49 deletions roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,6 +75,7 @@ using namespace RooStats ;
using namespace std ;

using namespace RooStats::HistFactory::Detail;
using namespace RooStats::HistFactory::Detail::MagicConstants;

ClassImp(RooStats::HistFactory::HistoToWorkspaceFactoryFast);

Expand Down Expand Up @@ -968,7 +969,7 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
RooArgList statFactorParams = ParamHistFunc::createParamSet(proto,
ParamSetPrefix,
theObservables,
0.0, 10.0);
defaultGammaMin, defaultStatErrorGammaMax);

ParamHistFunc statUncertFunc(statFuncName.c_str(), statFuncName.c_str(),
theObservables, statFactorParams );
Expand Down Expand Up @@ -1018,7 +1019,7 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
// We should change this to range from 0 to /inf
RooArgList shapeFactorParams = ParamHistFunc::createParamSet(proto,
funcParams,
theObservables, 0, 1000);
theObservables, defaultGammaMin, defaultShapeFactorGammaMax);

// Create the Function
ParamHistFunc shapeFactorFunc( funcName.c_str(), funcName.c_str(),
Expand Down Expand Up @@ -1062,7 +1063,7 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
cxcoutF(HistFactory) << "Cannot include Stat Error for histograms of more than 3 dimensions."
<< std::endl;
throw hf_exc();
} else {
}

// List of ShapeSys ParamHistFuncs
std::vector<string> ShapeSysNames;
Expand Down Expand Up @@ -1098,7 +1099,7 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
std::string funcParams = "gamma_" + shapeSys.GetName();
RooArgList shapeFactorParams = ParamHistFunc::createParamSet(proto,
funcParams,
theObservables, 0, 10);
theObservables, defaultGammaMin, defaultShapeSysGammaMax);

// Create the Function
ParamHistFunc shapeFactorFunc( funcName.c_str(), funcName.c_str(),
Expand All @@ -1125,8 +1126,7 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
systype = Constraint::Poisson;
}

double minShapeUncertainty = 0.0;
RooArgList shapeConstraints = createStatConstraintTerms(proto, constraintTermNames,
RooArgList shapeConstraints = createGammaConstraintTerms(proto, constraintTermNames,
*paramHist, shapeErrorHist,
systype,
minShapeUncertainty);
Expand All @@ -1141,8 +1141,6 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
sampleHistFuncs.push_back(proto.function(name));
}

} // End: NumObsVar == 1

} // End: !GetShapeSysList.empty()


Expand Down Expand Up @@ -1206,7 +1204,7 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
}

double statRelErrorThreshold = channel.GetStatErrorConfig().GetRelErrorThreshold();
RooArgList statConstraints = createStatConstraintTerms(proto, constraintTermNames,
RooArgList statConstraints = createGammaConstraintTerms(proto, constraintTermNames,
*chanStatUncertFunc, fracStatError.get(),
statConstraintType,
statRelErrorThreshold);
Expand Down Expand Up @@ -1754,9 +1752,8 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
}



RooArgList HistoToWorkspaceFactoryFast::
createStatConstraintTerms( RooWorkspace& proto, vector<string>& constraintTermNames,
createGammaConstraintTerms( RooWorkspace& proto, std::vector<std::string>& constraintTermNames,
ParamHistFunc& paramHist, const TH1* uncertHist,
Constraint::Type type, double minSigma ) {

Expand Down Expand Up @@ -1787,29 +1784,27 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
// Check that there are N elements
// in the RooArgList
if( numBins != numParams ) {
std::cout << "Error: In createStatConstraintTerms, encountered bad number of bins" << std::endl;
std::cout << "Error: In createGammaConstraintTerms, encountered bad number of bins" << std::endl;
std::cout << "Given histogram with " << numBins << " bins,"
<< " but require exactly " << numParams << std::endl;
throw hf_exc();
}

int TH1BinNumber = 0;
for( int i = 0; i < paramSet.getSize(); ++i) {

TH1BinNumber++;
// Get the relative sigmas from the hist
std::vector<double> relSigmas(paramSet.size());
for(int i = 0; i < paramSet.getSize(); ++i) {
relSigmas[i] = uncertHist->GetBinContent(i + 1);
}
configureConstrainedGammas(paramSet, relSigmas, minSigma);

while( uncertHist->IsBinUnderflow(TH1BinNumber) || uncertHist->IsBinOverflow(TH1BinNumber) ){
TH1BinNumber++;
}
for( int i = 0; i < paramSet.getSize(); ++i) {

RooRealVar& gamma = (RooRealVar&) (paramSet[i]);

cxcoutI(HistFactory) << "Creating constraint for: " << gamma.GetName()
<< ". Type of constraint: " << type << std::endl;

// Get the sigma from the hist
// (the relative uncertainty)
const double sigmaRel = uncertHist->GetBinContent(TH1BinNumber);
const double sigmaRel = relSigmas[i];

// If the sigma is <= 0,
// do cont create the term
Expand All @@ -1818,16 +1813,11 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
<< gamma.GetName()
<< " because sigma = " << sigmaRel
<< " (sigma<=0)"
<< " (TH1 bin number = " << TH1BinNumber << ")"
<< " (TH1 bin number = " << (i + 1) << ")"
<< std::endl;
gamma.setConstant(true);
continue;
}

// set reasonable ranges for gamma parameters
gamma.setMax( 1 + 5*sigmaRel );
gamma.setMin( 0. );

// Make Constraint Term
std::string constrName = std::string(gamma.GetName()) + "_constraint";
std::string nomName = std::string("nom_") + gamma.GetName();
Expand All @@ -1846,10 +1836,6 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo

// Make the constraint:
getOrCreate<RooGaussian>(proto, constrName, constrNom, gamma, constrSigma);

// Give reasonable starting point for pre-fit errors by setting it to the absolute sigma
// Mostly useful for pre-fit plotting.
gamma.setError(sigmaRel);
} else if( type == Constraint::Poisson ) {

double tau = 1/sigmaRel/sigmaRel; // this is correct Poisson equivalent to a Gaussian with mean 1 and stdev sigma
Expand All @@ -1870,30 +1856,13 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
// Type 2 : RooPoisson
getOrCreate<RooPoisson>(proto, constrName, constrNom, constrMean).setNoRounding(true);

if (std::string(gamma.GetName()).find("gamma_stat") != std::string::npos) {
// Give reasonable starting point for pre-fit errors.
// Mostly useful for pre-fit plotting.
gamma.setError(sigmaRel);
}

} else {

std::cout << "Error: Did not recognize Stat Error constraint term type: "
<< type << " for : " << paramHist.GetName() << std::endl;
throw hf_exc();
}

// If the sigma value is less
// than a supplied threshold,
// set the variable to constant
if( sigmaRel < minSigma ) {
cxcoutW(HistFactory) << "Warning: Bin " << i << " = " << sigmaRel
<< " and is < " << minSigma
<< ". Setting: " << gamma.GetName() << " to constant"
<< std::endl;
gamma.setConstant(true);
}

constraintTermNames.push_back( constrName );
ConstraintTerms.add( *proto.pdf(constrName) );

Expand Down
27 changes: 8 additions & 19 deletions roofit/hs3/src/JSONFactories_HistFactory.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,9 @@

using RooFit::Detail::JSONNode;

using namespace RooStats::HistFactory::Detail;
using namespace RooStats::HistFactory::Detail::MagicConstants;

namespace {

inline void writeAxis(JSONNode &axis, RooRealVar const &obs)
Expand Down Expand Up @@ -141,8 +144,6 @@ std::string toString(TClass *c)
return "unknown";
}

using namespace RooStats::HistFactory::Detail;

RooRealVar &createNominal(RooWorkspace &ws, std::string const &parname, double val, double min, double max)
{
RooRealVar &nom = getOrCreate<RooRealVar>(ws, "nom_" + parname, val, min, max);
Expand All @@ -165,20 +166,6 @@ RooAbsPdf &getConstraint(RooWorkspace &ws, const std::string &sysname, const std
return getOrCreate<RooGaussian>(ws, constraintName(sysname), *constrParam, *ws.var("nom_" + pname), 1.);
}

void setMCStatGammaRanges(RooArgList const &gammas, std::vector<double> const &errs, double statErrThresh)
{
// Set a reasonable range for gamma and set constant NPs which are below the
// MC stat threshold, and remove them from the np list.
for (size_t i = 0; i < errs.size(); ++i) {
auto g = static_cast<RooRealVar *>(gammas.at(i));
g->setMax(1 + 5 * errs[i]);
g->setError(errs[i]);
if (errs[i] < statErrThresh) {
g->setConstant(true); // all negative errs are set constant
}
}
}

ParamHistFunc &createPHF(const std::string &sysname, const std::string &phfname, const std::vector<double> &vals,
RooJSONFactoryWSTool &tool, RooArgList &constraints, const RooArgSet &observables,
const std::string &constraintType, RooArgList &gammas, double gamma_min, double gamma_max)
Expand Down Expand Up @@ -319,7 +306,9 @@ bool importHistSample(RooJSONFactoryWSTool &tool, RooDataHist &dh, RooArgSet con
}
RooArgList gammas;
std::string constraint(mod["constraint"].val());
shapeElems.add(createPHF(sysname, funcName, vals, tool, constraints, varlist, constraint, gammas, 0, 1000));
shapeElems.add(createPHF(sysname, funcName, vals, tool, constraints, varlist, constraint, gammas,
defaultGammaMin, defaultShapeSysGammaMax));
configureConstrainedGammas(gammas, vals, minShapeUncertainty);
} else if (modtype == "custom") {
RooAbsReal *obj = ws.function(sysname);
if (!obj) {
Expand Down Expand Up @@ -423,8 +412,8 @@ class HistFactoryImporter : public RooFit::JSONIO::Importer {

RooArgList gammas;
mcStatObject = &createPHF("stat_" + phfName, "mc_stat_" + phfName, vals, *tool, constraints, observables,
statErrType, gammas, 0, 10);
setMCStatGammaRanges(gammas, errs, statErrThresh);
statErrType, gammas, defaultGammaMin, defaultStatErrorGammaMax);
configureConstrainedGammas(gammas, errs, statErrThresh);
}

int idx = 0;
Expand Down

0 comments on commit 3c68044

Please sign in to comment.