#include <iostream>
#include <fstream>
#include <sstream>
#include <cstdlib>
#include <vector>
#include <iterator>
#include <algorithm>
#include "SFramework/TSModelFactory.h"
#include "RooStats/HistFactory/PiecewiseInterpolation.h"
#include "RooStats/HistFactory/FlexibleInterpVar.h"
#include "RooStats/ModelConfig.h"
#include "RooRealVar.h"
#include "RooRealSumPdf.h"
#include "QFramework/TQStringUtils.h"
#include "QFramework/TQTaggable.h"
#include "QFramework/TQIterator.h"
#include "QFramework/TQPathManager.h"
#include "QFramework/TQUtils.h"
#include "TSystem.h"
#include "TMath.h"
#include "TFile.h"
#ifdef __cpp_generic_lambdas
#if __cpp_generic_lambdas >= 201304
#define HAS_GENERIC_LAMBDAS
#endif
#endif
#include "SFramework/TSIterator.h"
using namespace RooStats::HistFactory;
ClassImp(TSModelFactory)
#ifdef HAS_GENERIC_LAMBDAS
namespace {
namespace SFINAE {
class true_val {char dummy[1];};
class false_val {char dummy[2];};
template <typename T>
class hasStackLabel {
public:
template <typename C> static true_val test(decltype(&C::SetStackLabel));
template <typename C> static false_val test(...);
enum { value = sizeof(test<T>(0)) == sizeof(true_val) };
};
template <typename T, typename F> auto static_if(std::true_type, T t, F ) { return t; }
template <typename T, typename F> auto static_if(std::false_type, T , F f) { return f; }
template <bool B, typename T, typename F> auto static_if(T t, F f) { return static_if(std::integral_constant<bool, B>{}, t, f); }
template <bool B, typename T> auto static_if(T t) { return static_if(std::integral_constant<bool, B>{}, t, [](auto&&...){}); }
template <class C> inline void trySetStackLabel(C& staterr, const char* stackLabel){
static_if<hasStackLabel<C>::value>
([&](auto& s) {
s.SetStackLabel(stackLabel);
})
(staterr);
}
}
}
#endif
TSModelFactory::TSModelFactory() {
fShowInfo = false;
fShowError = true;
fShowWarn = false;
}
TSModelFactory::~TSModelFactory() {
}
void TSModelFactory::info(TString message) {
if (fShowInfo) {
std::cout << "SFramework/TSModelFactory: " << message.Data() << std::endl;
}
}
void TSModelFactory::error(TString message) {
if (fShowError) {
std::cout << "SFramework/TSModelFactory: " << TQStringUtils::makeBoldRed(TString("ERROR: ") + message) << std::endl;
}
}
void TSModelFactory::warn(TString message) {
if (fShowWarn) {
std::cout << "SFramework/TSModelFactory: " << TQStringUtils::makeBoldYellow(TString("WARNING: ") + message) << std::endl;
}
}
Bool_t TSModelFactory::createXML(TQFolder * model, const TString& xmlDir) {
RooStats::HistFactory::Measurement * ms = createMeasurement(model);
if (ms) {
Bool_t success = createXML(ms, xmlDir, model->getTagStringDefault(".HistogramFile", TQFolder::concatPaths(xmlDir,"histograms.root")));
delete ms;
return success;
} else {
return false;
}
}
Bool_t TSModelFactory::writeHistograms(TQFolder* model, const TString& hfname){
TQFolder * histos = model->getFolder(".Histograms");
if (!histos) return false;
TString histosFilename = model->replaceInText(hfname);
histosFilename = TQPathManager::getPathManager()->getTargetPath(histosFilename);
if (!TQUtils::ensureDirectoryForFile(histosFilename)) return false;
TFile * histosFile = TFile::Open(histosFilename, "RECREATE", histosFilename.Data());
if (!histosFile || !histosFile->IsOpen()) {
if (histosFile) {
delete histosFile;
}
return false;
}
TString targetName = TString::Format("%s_histograms", model->GetName());
model->setTagString(".HistogramFile", histosFilename);
model->setTagString(".HistogramPathPrefix", targetName);
TDirectory* target = histosFile->mkdir(targetName);
TQFolder* nominal = histos->getFolder("Nominal");
if(nominal){
nominal->writeDirectory(target);
} else {
return false;
}
TQFolder* shapesysts = histos->getFolder("Systematics");
if(shapesysts){
shapesysts->writeDirectory(target);
}
histosFile->Write();
histosFile->Close();
delete histosFile;
return true;
}
Bool_t TSModelFactory::createXML(RooStats::HistFactory::Measurement * measurement, const TString& xmlDir, const TString& histosFile) {
if (measurement) {
measurement->PrintXML(xmlDir.Data());
char* rootsys = std::getenv("ROOTSYS");
TString s = "cp ";
s.Append(TQFolder::concatPaths(rootsys,"etc","HistFactorySchema.dtd"));
s.Append(" ");
s.Append(xmlDir);
s.Append(+"/");
info(s);
TList* l = TQUtils::execute(s);
delete l;
l = TQUtils::execute ("sed -i -e 's|<Data[ ]*HistoName=\"\"[ ]*InputFile=\"\"[ ]*HistoPath=\"\"[ ]*/>||' "+xmlDir+"/*.xml");
delete l;
l = TQUtils::execute ("sed -i -e 's|InputFile=\"[^\"]*\"|InputFile=\""+histosFile+"\"|' "+xmlDir+"/*.xml");
delete l;
TString cwd = gSystem->pwd();
TQStringUtils::ensureTrailingText(cwd,"/");
l = TQUtils::execute ("sed -i -e 's|<Input>./"+cwd+"|<Input>|' "+xmlDir+"/*.xml");
delete l;
return true;
} else {
return false;
}
}
RooStats::HistFactory::Measurement * TSModelFactory::createMeasurement(TQFolder * model) {
if (!model) {
error("createMeasurement(): No model definition provided");
return 0;
}
fModel = model;
info("createMeasurement(...): Creating Measurement object...");
RooStats::HistFactory::Measurement * measurement = new RooStats::HistFactory::Measurement(model->GetName(), model->getTagStringDefault("Title", model->GetTitle()));
TString outputFilePrefix = model->getTagStringDefault("OutputFilePrefix", model->GetName());
measurement->SetOutputFilePrefix(outputFilePrefix.Data());
Bool_t exportOnly = true;
if (model->getTagBool("ExportOnly", exportOnly))
measurement->SetExportOnly(exportOnly);
auto pois = model->getTagVString("POI");
std::sort(pois.begin(),pois.end());
for(auto poi:pois){
measurement->AddPOI(poi.Data());
}
Double_t lumi = 1.;
if (model->getTagDouble("Lumi", lumi))
measurement->SetLumi(lumi);
Double_t lumiRelErr = 0.;
if (model->getTagDouble("LumiRelErr", lumiRelErr))
measurement->SetLumiRelErr(lumiRelErr);
TQIterator itrCh(model->getListOfFolders("?"), "!.*", true);
while (itrCh.hasNext()) {
addChannel((TQFolder*)itrCh.readNext(), measurement);
}
TQIterator itrParam(model->getListOfFolders("ParamSettings/?"), "!.*", true);
while (itrParam.hasNext()) {
TQFolder * folder = (TQFolder*)itrParam.readNext();
TQTaggable tagsParam(folder);
Bool_t isConst = false;
if (tagsParam.getTagBool("Const", isConst) && isConst) {
measurement->AddConstantParam(folder->GetName());
}
if (tagsParam.getTagBoolDefault("POI", false)) {
measurement->AddPOI(folder->GetName());
}
TString constraint;
if (tagsParam.getTagString("ConstraintTerm", constraint)) {
if (constraint.CompareTo("Uniform", TString::kIgnoreCase) == 0) {
measurement->AddUniformSyst(folder->GetName());
} else if (constraint.CompareTo("Gamma", TString::kIgnoreCase) == 0) {
measurement->AddGammaSyst(folder->GetName(), 1.);
} else if (constraint.CompareTo("LogNorm", TString::kIgnoreCase) == 0) {
measurement->AddLogNormSyst(folder->GetName(), 1.);
} else if (constraint.CompareTo("Gaussian", TString::kIgnoreCase) != 0) {
warn(TString::Format("createMeasurement(): Unknown constraint term '%s' "
"for parameter '%s'", constraint.Data(), folder->GetName()));
}
}
}
fModel = 0;
return measurement;
}
Bool_t TSModelFactory::reviseWorkspace(RooWorkspace * workspace, TQFolder * model) {
if (!workspace) {
error("reviseWorkspace(): No valid workspace provided");
return false;
}
if (!model) {
error("reviseWorkspace(): No model definition provided");
return false;
}
RooStats::ModelConfig* mc = dynamic_cast<RooStats::ModelConfig*>(workspace->obj("ModelConfig"));
if (!mc) {
error("reviseWorkspace(): No ModelConfig available");
return false;
}
workspace->SetTitle(model->getTagStringDefault("Title", model->GetTitle()));
TQFolder* paramSettings = model->getFolder("ParamSettings");
RooArgSet thenps;
const RooAbsCollection* nps = mc->GetNuisanceParameters();
if(!nps || nps->getSize() < 1){
warn("no nuisance parameters set, setting manually");
RooArgList allVars(workspace->allVars());
RooRealVarIterator itr(allVars.createIterator());
while (itr.hasNext()) {
RooRealVar * var = itr.readNext();
if(!var) continue;
bool isNP = TQStringUtils::matches(var->GetName(),"alpha_*") || TQStringUtils::matches(var->GetName(),"gamma_stat_*");
if(paramSettings){
TQFolder* p = paramSettings->getFolder(var->GetName());
if(p){
p->getTagBool("isNP",isNP);
}
}
if(isNP){
thenps.add(*var);
}
}
} else {
thenps.add(*nps);
}
auto pois = model->getTagVString("POI");
RooArgSet thepois;
if(mc->GetParametersOfInterest()){
thepois.add(*mc->GetParametersOfInterest());
}
for(auto poi:pois){
for(const auto& var:workspace->allVars()){
if(TQStringUtils::matches(var->GetName(),poi)){
thepois.add(*var);
}
}
}
RooArgSet finalnps;
for(const auto& np:thenps){
if(!thepois.find(np->GetName())) finalnps.add(*np);
}
mc->SetParametersOfInterest(thepois);
mc->SetNuisanceParameters(finalnps);
if(paramSettings){
TQIterator itrParam(paramSettings->getListOfFolders("?"), "!.*", true);
while (itrParam.hasNext()) {
TQFolder * folder = (TQFolder*)itrParam.readNext();
TString name = folder->GetName();
RooRealVar * var = workspace->var(name.Data());
if (!var) {
warn(TString::Format("reviseWorkspace(): Failed to find variable '%s' "
"in workspace. Skipping ...", name.Data()));
continue;
}
TQTaggable tagsParam(folder);
TString title;
if (tagsParam.getTagString("Title", title)) {
var->SetTitle(title.Data());
}
Int_t interpCode = 4;
if (tagsParam.getTagInteger("InterpCode", interpCode)) {
TQIterator itr(var->clientIterator(), true);
while (itr.hasNext()) {
TObject * obj = itr.readNext();
if (obj->InheritsFrom(FlexibleInterpVar::Class())) {
((FlexibleInterpVar*)obj)->setInterpCode(*var, interpCode);
} else if (obj->InheritsFrom(PiecewiseInterpolation::Class())) {
((PiecewiseInterpolation*)obj)->setInterpCode(*var, interpCode);
} else {
error(TString::Format("reviseWorkspace(): Failed to set interpolation code "
"on variable '%s'. Skipping ...", var->GetName()));
}
}
}
TString desc;
if(tagsParam.getTagString("Description",desc)){
var->setAttribute(desc.Data(),true);
}
std::vector<TString> attributes = tagsParam.getTagVString("Attributes");
for(const auto& attr:attributes){
var->setAttribute(attr.Data());
}
}
}
if(model->getTagBoolDefault("isBinned",true)){
RooFIter iter = workspace->components().fwdIterator() ;
RooAbsArg* arg = NULL;
while((arg = iter.next())) {
if (arg->IsA() == RooRealSumPdf::Class()) {
arg->setAttribute("BinnedLikelihood");
}
}
}
RooRealVar* lumi = workspace->var("Lumi");
if(lumi){
lumi->setConstant(true);
}
TList * lConfig = model->exportToText(false);
workspace->import(*lConfig, ".Model", true);
return true;
}
Bool_t TSModelFactory::addChannel(TQFolder * definition, RooStats::HistFactory::Measurement * measurement) {
TString channelName = definition->GetName();
if (!TQStringUtils::removeLeadingText(channelName, "Channel."))
return false;
TString histosfile = definition->getTagStringDefault("~.HistogramFile", "");
TString histospath = TQFolder::concatPaths(definition->getTagStringDefault("~.HistogramPathPrefix", ""),"Nominal",channelName.Data());
RooStats::HistFactory::Channel * channel = new Channel(channelName.Data(),histospath.Data());
channel->SetStatErrorConfig(
definition->getTagDoubleDefault("StatRelErrorThreshold", 0.05),
definition->getTagStringDefault("StatConstraintType", "Poisson").Data());
TQIterator itr(definition->getListOfFolders("?"), "!.*", true);
while (itr.hasNext()) {
TQFolder * folder = (TQFolder*)itr.readNext();
if (addSample(folder, channel)) {
info(TString::Format("Added as Sample '%s'", folder->GetName()));
}
}
measurement->AddChannel(*channel);
return true;
}
Bool_t TSModelFactory::addSample(TQFolder * definition, RooStats::HistFactory::Channel * channel) {
TString sampleName = definition->GetName();
if (!TQStringUtils::removeLeadingText(sampleName, "Sample.")) {
return false;
}
TString sampleType = definition->getTagStringDefault("Type","B");
TString histogramLocation;
if(!definition->getTagString("Histo",histogramLocation)){
error(TString::Format("unable to find histogram for '%s'",definition->getPath().Data()));
return false;
}
TString histogramFile = definition->getTagStringDefault("~.HistogramFile","histFactor_tmp.root");
TString histogramPathPrefix = definition->getTagStringDefault("~.HistogramPathPrefix","");
TString histoPath = TQFolder::concatPaths(histogramPathPrefix, histogramLocation);
TString histoName = TQFolder::getPathTail(histoPath);
TQStringUtils::ensureTrailingText(histoPath,"/");
TString histoFilePath = TQPathManager::getPathManager()->getTargetPath(histogramFile).c_str();
if (sampleType != "D"){
Sample * sample = new Sample(sampleName.Data(),histoName.Data(), histoFilePath.Data(), histoPath.Data());
sample->SetNormalizeByTheory(definition->getTagBoolDefault("NormalizeByTheory", true));
if (definition->getTagBoolDefault("ActivateStatError", false)) {
sample->GetStatError().Activate(true);
#ifdef HAS_GENERIC_LAMBDAS
TString stackLabel;
if(definition->getTagString("StatErrorStackLabel",stackLabel)){
::SFINAE::trySetStackLabel(sample->GetStatError(),stackLabel.Data());
}
#endif
if (definition->hasTag("StatErrorHisto")) {
TString statErrorHistoFile;
TString statErrorHistoName;
TString statErrorHistoPath;
splitHistoLocation(definition->getTagStringDefault("StatErrorHisto"),statErrorHistoFile, statErrorHistoName, statErrorHistoPath);
TString statErrorHistoFilePath = TQPathManager::getPathManager()->getTargetPath(statErrorHistoFile).c_str();
sample->GetStatError().SetUseHisto(true);
sample->GetStatError().SetHistoName(statErrorHistoName.Data());
sample->GetStatError().SetHistoName(statErrorHistoFilePath.Data());
sample->GetStatError().SetHistoName(statErrorHistoPath.Data());
}
} else {
sample->GetStatError().Activate(false);
}
TQIterator itr(definition->getListOfFolders("?"), "!.*", true);
while (itr.hasNext()) {
TQFolder * folder = (TQFolder*)itr.readNext();
TString name = folder->GetName();
if (addNormFactor(folder, sample)) {
info(TString::Format("Added as NormFactor '%s'", folder->GetName()));
} else if (addOverallSys(folder, sample)) {
info(TString::Format("Added as OverallSys '%s'", folder->GetName()));
} else if (addHistoSys(folder, sample)) {
info(TString::Format("Added as HistoSys '%s'", folder->GetName()));
} else if (addShapeSys(folder, sample)) {
info(TString::Format("Added as ShapeSys '%s'", folder->GetName()));
} else if (addShapeFactor(folder, sample)) {
info(TString::Format("Added as ShapeFactor '%s'", folder->GetName()));
} else if (addHistoFactor(folder, sample)) {
info(TString::Format("Added as HistFactor '%s'", folder->GetName()));
} else {
error(TString::Format("Unrecognized element '%s' "
"in sample '%s'", name.Data(), sampleName.Data()));
}
}
channel->AddSample(*sample);
} else {
TString dataName = definition->getTagStringDefault("dataset", sampleName == "Data" ? "obsData" : sampleName.Data());
if(dataName == "obsData"){
channel->SetData(histoName.Data(), histoFilePath.Data(), histoPath.Data());
} else {
auto data = RooStats::HistFactory::Data(histoName.Data(), histoFilePath.Data(), histoPath.Data());
data.SetName(dataName.Data());
channel->AddAdditionalData(data);
}
}
return true;
}
Bool_t TSModelFactory::addNormFactor(TQFolder * definition, Sample * sample) {
TString normFactorName = definition->GetName();
if (!TQStringUtils::removeLeadingText(normFactorName, "NormFactor.")) {
return false;
}
sample->AddNormFactor(normFactorName.Data(),
definition->getTagDoubleDefault("Val", 1.),
definition->getTagDoubleDefault("Low", 1.),
definition->getTagDoubleDefault("High", 1.));
TQFolder * paramSettings = fModel->getFolder(TQFolder::concatPaths("ParamSettings", normFactorName) + "+");
Bool_t isConst = false;
if (definition->getTagBool("Const", isConst)) {
TQTaggable tagsParam(paramSettings);
if (tagsParam.getTagBoolDefault("Const", isConst) == isConst) {
paramSettings->setTagBool("Const", isConst);
} else {
paramSettings->setTagBool(".Const.Recessive", isConst);
}
}
if(definition->getTagBoolDefault("POI",false)){
paramSettings->setTagBool("POI",true);
}
return true;
}
Bool_t TSModelFactory::addShapeFactor(TQFolder * definition, Sample * sample) {
TString shapeFactorName = definition->GetName();
if (!TQStringUtils::removeLeadingText(shapeFactorName, "ShapeFactor.")) {
return false;
}
sample->AddShapeFactor(shapeFactorName.Data());
TQFolder * paramSettings = fModel->getFolder(TQFolder::concatPaths("ParamSettings", shapeFactorName) + "+");
if(definition->getTagBoolDefault("POI",false)){
paramSettings->setTagBool("POI",true);
}
return true;
}
Bool_t TSModelFactory::addHistoFactor(TQFolder * definition, Sample * sample) {
TString histoFactorName = definition->GetName();
if (!TQStringUtils::removeLeadingText(histoFactorName, "HistoFactor.")) {
return false;
}
TString histogramLocationHigh,histogramLocationLow;
if(!definition->getTagString("HistoHigh",histogramLocationHigh) || !definition->getTagString("HistoLow",histogramLocationLow)){
error(TString::Format("no histogram assigned to '%s'",definition->getPath().Data()));
return false;
}
TString histogramFile = definition->getTagStringDefault("~.HistogramFile","histFactor_tmp.root");
TString histogramPathPrefix = definition->getTagStringDefault("~.HistogramPathPrefix","");
TString histoPathLow = TQFolder::concatPaths(histogramPathPrefix, histogramLocationLow);
TString histoPathHigh = TQFolder::concatPaths(histogramPathPrefix, histogramLocationHigh);
TString histoNameLow = TQFolder::getPathTail(histoPathLow);
TString histoNameHigh = TQFolder::getPathTail(histoPathHigh);
TQStringUtils::ensureTrailingText(histoPathLow,"/");
TQStringUtils::ensureTrailingText(histoPathHigh,"/");
TString histoFilePath = TQPathManager::getPathManager()->getTargetPath(histogramFile).c_str();
sample->AddHistoFactor(histoFactorName.Data(),
histoNameLow.Data(), histoFilePath.Data(), histoPathLow.Data(),
histoNameHigh.Data(), histoFilePath.Data(), histoPathHigh.Data());
return true;
}
Bool_t TSModelFactory::addOverallSys(TQFolder * definition, Sample * sample) {
TString overallSysName = definition->GetName();
if (!TQStringUtils::removeLeadingText(overallSysName, "OverallSys.")) {
return false;
}
sample->AddOverallSys(overallSysName.Data(),
definition->getTagDoubleDefault("Low", 1.),
definition->getTagDoubleDefault("High", 1.));
return true;
}
Bool_t TSModelFactory::addHistoSys(TQFolder * definition, Sample * sample) {
TString histoSysName = definition->GetName();
if (!TQStringUtils::removeLeadingText(histoSysName, "HistoSys.")) {
return false;
}
TString histogramLocationHigh,histogramLocationLow;
if(!definition->getTagString("HistoHigh",histogramLocationHigh) || !definition->getTagString("HistoLow",histogramLocationLow)){
error(TString::Format("no histogram assigned to '%s'",definition->getPath().Data()));
return false;
}
TString histogramFile = definition->getTagStringDefault("~.HistogramFile","histFactor_tmp.root");
TString histogramPathPrefix = definition->getTagStringDefault("~.HistogramPathPrefix","");
TString histoPathLow = TQFolder::concatPaths(histogramPathPrefix, histogramLocationLow);
TString histoPathHigh = TQFolder::concatPaths(histogramPathPrefix, histogramLocationHigh);
TString histoNameLow = TQFolder::getPathTail(histoPathLow);
TString histoNameHigh = TQFolder::getPathTail(histoPathHigh);
TQStringUtils::ensureTrailingText(histoPathLow,"/");
TQStringUtils::ensureTrailingText(histoPathHigh,"/");
TString histoFilePath = TQPathManager::getPathManager()->getTargetPath(histogramFile).c_str();
sample->AddHistoSys(histoSysName.Data(),
histoNameLow.Data(), histoFilePath.Data(), histoPathLow.Data(),
histoNameHigh.Data(), histoFilePath.Data(), histoPathHigh.Data());
return true;
}
Bool_t TSModelFactory::addShapeSys(TQFolder * definition, Sample * sample) {
TString shapeSysName = definition->GetName();
if (!TQStringUtils::removeLeadingText(shapeSysName, "ShapeSys.")) {
return false;
}
TString histogramLocation;
if(!definition->getTagString("Histo",histogramLocation)){
error(TString::Format("unable to find histogram for '%s'",definition->getPath().Data()));
return false;
}
TString histogramFile;
TString histogramPathPrefix;
definition->getTagString("~.HistogramFile", histogramFile, true);
definition->getTagString("~.HistogramPathPrefix", histogramPathPrefix, true);
bool usePoisson = definition->getTagBoolDefault("usePoisson",false);
histogramLocation = histogramFile + ":"
+ TQFolder::concatPaths(histogramPathPrefix, histogramLocation);
TString histoFile;
TString histoName;
TString histoPath;
splitHistoLocation(histogramLocation, histoFile, histoName, histoPath);
sample->AddShapeSys(shapeSysName.Data(), ( usePoisson ? RooStats::HistFactory::Constraint::Type::Poisson : RooStats::HistFactory::Constraint::Type::Gaussian ) , histoName.Data(), TQPathManager::getPathManager()->getTargetPath(histoFile).c_str(), histoPath.Data());
return true;
}
void TSModelFactory::splitHistoLocation(TString input,
TString &histoFile, TString &histoName, TString &histoPath) {
histoPath = input;
histoFile = TQStringUtils::readPrefix(histoPath, ":");
histoName = TQFolder::getPathTail(histoPath);
if (!histoPath.IsNull())
histoPath.Append("/");
}