-
Carsten Burgard authoredCarsten Burgard authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
ActionScaleSystematics.cxx 10.61 KiB
#include "QFramework/TQFolder.h"
#include "QFramework/TQIterator.h"
#include "QFramework/TQUtils.h"
#include "QFramework/TQStringUtils.h"
#include "QFramework/TQHistogramUtils.h"
#include "SFramework/TSStatisticsManager.h"
#include "SFramework/TSUtils.h"
#include "SFramework/TSHistogramExporter.h"
#include "RooStats/ModelConfig.h"
#include "RooCategory.h"
#include "RooBinning.h"
#include "RooDataSet.h"
#include "RooRealSumPdf.h"
#include "RooProdPdf.h"
#include "RooArgSet.h"
#include "RooProduct.h"
#include "RooRealVar.h"
#include "RooAbsReal.h"
#include "RooFitResult.h"
#include "RooAddition.h"
#include "RooSimultaneous.h"
#include "RooPoisson.h"
#include "RooConstVar.h"
#include "RooHistFunc.h"
#include "RooDataHist.h"
#include "RooStats/HistFactory/PiecewiseInterpolation.h"
#include "RooStats/HistFactory/FlexibleInterpVar.h"
#include "RooStats/HistFactory/ParamHistFunc.h"
#include "RooAbsDataStore.h"
#include "TFile.h"
#include "RooCmdConfig.h"
#include <map>
#include <utility>
/*<cafdoc name=ScaleSystematics>
ScaleSystematics
===========================
This action scales the effect of certain parameters (systematic uncertainties)
in a workspace by an arbitrary factor.
The Action scales the effect of all (matching) components of
`PiecewiseInterpolation` and `RooStats::HistFactory::FlexibleInterpVar`
instances in a workspace/pdf.
Note: This has been designed to work with workspaces like the HWW one when
created using `SFramework`/`HistFactory`.
Usage:
---------------
```
+ScaleSystematics{
+NameOfWorkspace {
<pdf = "nameOfPdfToEdit"> #defaults to "simPdf"
+someParameter{
<name = "expression*ToMatch", scale = 0.5, invert=false, sqrt=false>
#scale the effects of parameter(s) with names matching the expression given as
#'name' by the given scaling factor. If 'invert' is true, yields
#are scaled by 1/'scale' instead (default is false).
#if 'sqrt' is true the sqrt of the scaling factor is used (or its
#inverse)
}
+otherSyst {
<scale = 3, invert=true>
#if 'name' is not specified, the name of the folder is used (here:
# "otherProcess")
}
}
}
```
The process name is matched against the name of the relevant `RooAbsArg`
instances. If their names are not known, try printing a workspace in
cling (`root` shell): `ws->Print("t")`
</cafdoc> */
namespace TSBaseActions {
class ScaleSystematics : public TSStatisticsManager::Action {
bool doApply(const char* name, const TString& filter, const TString& veto) const {
return TQStringUtils::matchesFilter(name, filter,",") && !TQStringUtils::matchesFilter(name,veto,",");
}
void processPiecewiseInterpolation(const PiecewiseInterpolation* pw, const TString& filter, const TString& veto, const double scale) const {
RooAbsArg* nominalAbs = TSUtils::getNominalMember(pw);
if (!nominalAbs) return;
RooListProxy* paramSet = TSUtils::getParamSetMember(pw);
if (!paramSet) return;
RooListProxy* highSet = TSUtils::getHighSetMember(pw);
if (!highSet) return;
RooListProxy* lowSet = TSUtils::getLowSetMember(pw);
if (!lowSet) return;
auto paramIter(TSUtils::toVector<RooAbsReal*>(*paramSet));
auto highIter(TSUtils::toVector<RooAbsReal*>(*highSet)) ;
auto lowIter(TSUtils::toVector<RooAbsReal*>(*lowSet)) ;
RooHistFunc* nominalHist = dynamic_cast<RooHistFunc*>(nominalAbs);
for(size_t i=0; i<paramIter.size(); ++i){
auto parAbs = paramIter[i];
auto lowAbs = lowIter[i];
auto highAbs = highIter[i];
if (!doApply(parAbs->GetName(), filter, veto)) continue; //only process if this parameter matches our filter expression
std::cout<<"Scaling PiecewiseInterpolation component '"<<parAbs->GetName()<<"'"<<std::endl;
//assumption: this construct only appears in use with RooHistFunc s, other cases would still need to be implemented
if (nominalHist) {
RooHistFunc* highHist = dynamic_cast<RooHistFunc*>(highAbs);
RooHistFunc* lowHist = dynamic_cast<RooHistFunc*>(lowAbs);
if (!highHist || !lowHist) { //this should never happen as it would mean that nominal and variations are different types of objects....
throw std::runtime_error(TString::Format("Central value (nominal) '%s' is a RooHistFunc but high and/or low variations ('%s', '%s') are not... It seems your workspace has some severe issues!",nominalHist->GetName(), highAbs->GetName(), lowAbs->GetName()).Data());
}
TSUtils::scaleDifferenceToNominal(TSUtils::getDataHist(highHist),TSUtils::getDataHist(nominalHist),scale);
TSUtils::scaleDifferenceToNominal(TSUtils::getDataHist(lowHist),TSUtils::getDataHist(nominalHist),scale);
} else { //to support other data types, add treatment here (in addition to attempted dynamic_casts above):
std::cout<<"Encountered unsuported type of PiecewiseInterpolation component: '"<<nominalAbs->GetName()<<"' is of type '"<<nominalAbs->IsA()->GetName()<<"' instance"<<std::endl;
}
}
}
void processFlexibleInterpVar(const RooStats::HistFactory::FlexibleInterpVar* interp, const TString& filter, const TString& veto, const double scale) const {
if (!interp) return;
Double_t nominal = TSUtils::getNominalMember(interp);
std::vector<double>* high = TSUtils::getHighMember(interp);
std::vector<double>* low = TSUtils::getLowMember(interp);
RooListProxy* paramList = TSUtils::getParamListMember(interp);
if (!high || !low || !paramList) {
throw std::runtime_error(TString::Format("Some members required for processing of FlexibleInterpVar '%s' could not be retrieved, please check your workspace/pdf!",interp->GetName()).Data());
return;
}
if (high->size() != low->size()) {
throw std::runtime_error(TString::Format("Low and High variation vectors of FlexibleInterpVar '%s' are of different length, please check your workspace/pdf!",interp->GetName()).Data());
return;
}
int i=0;
ROOFIT_ITERATE(*paramList,RooAbsArg,param){
//ignore this high/low combination if the parameter does not match our selection criteria
if (!doApply(param->GetName(), filter, veto)) continue;
std::cout<<"Scaling effect of FlexibleInterpVar parameter '"<<param->GetName()<<"' by factor "<<scale<<std::endl;
//apply the actual scaling
high->at(i) = nominal + scale * (high->at(i)-nominal);
low->at(i) = nominal + scale * (low->at(i)-nominal);
++i;
}
}
bool execute(TQFolder * config) const override {
if(!config){
throw std::runtime_error("received NULL pointer as config object!");
}
//retrieve workspace
RooWorkspace * ws = dynamic_cast<RooWorkspace*>(workspaces()->getObject(config->GetName()));
if(!ws){
manager->error(TString::Format("unable to load workspace '%s': no such workspace!",config->GetName()));
return false;
}
//retrieve pdf to edit
TString pdfName = config->getTagStringDefault("pdfName","simPdf");
RooAbsPdf* pdf = ws->pdf(pdfName);
if (!pdf){
manager->error(TString::Format("unable to retrieve pdf '%s' from workspace '%s': no such pdf!",pdfName.Data(),config->GetName()));
return false;
}
//this is a "manual" hack to get all nominal RooHistFuncs
//std::set<const RooAbsArg*>* allHistFuncs = TSUtils::getContributionsRecursive(simpdf, RooHistFunc::Class(), false);
std::set<const RooAbsArg*>* allPiecewise = TSUtils::getContributionsRecursive(pdf, PiecewiseInterpolation::Class(), false);
for (const auto* pw : *allPiecewise) {
if (!pw) continue;
pw->Print();
}
std::set<const RooAbsArg*>* allFlexInterp = TSUtils::getContributionsRecursive(pdf, RooStats::HistFactory::FlexibleInterpVar::Class(), false);
for (const auto* pw : *allFlexInterp) {
if (!pw) continue;
pw->Print();
}
TQFolderIterator operationsItr(config->getListOfFolders("?"),true); //retrieve the individual scaling operations we should perform
while (operationsItr.hasNext()) {
TQFolder* op = operationsItr.readNext();
//read some basic options for this part of the payload
//@tag: [name] Possibly comma separated list of (wildcarded) expressions the name of a systematic (nuissance parameter) needs to match for its effects to be modified. Defaults to name of the folder it is read from.
TString pattern = op->getTagStringDefault("name",op->GetName());
//@tag: [exclude] If the name of a systematic (nuissance parameter) matches this (possibly comma separated list of) expression(s) it is not modified even if it matches the 'name' pattern(s). Default: ''
TString vetoPattern = op->getTagStringDefault("exclude","");
//@tag: [scale] Value by which the effect of the selected systematics should be changed. Default: 1.0
double scale = op->getTagDoubleDefault("scale",1.);
//@tag: [invert] If set to true will multiply effects by '1/scale' instead of multiplying them by 'scale'
if (scale>0 && op->getTagBoolDefault("invert",false)) {
scale = 1./scale;
std::cout<<"inverting scale factor to "<<scale<<std::endl;
}
//@tag: [sqrt] If set the effect of selected parameters is multiplied by 'sqrt(scale)' instead of 'scale' (can be used together with 'invert')
if (op->getTagBoolDefault("sqrt",false)) {
scale = sqrt(scale);
std::cout<<"will scale by sqrt(scale) = "<<scale<<std::endl;
}
for (const RooAbsArg* pw_ : *allPiecewise) {
if (!pw_) continue;
const PiecewiseInterpolation* pw = dynamic_cast<const PiecewiseInterpolation*>(pw_);
if (!pw) {
std::cout<<"Skipping object '"<<pw_->GetName()<<"': not a PiecewiseInterpolation"<<std::endl;
continue;
}
processPiecewiseInterpolation(pw, pattern, vetoPattern, scale);
}
for (const RooAbsArg* pw_ : *allFlexInterp) {
if (!pw_) continue;
const RooStats::HistFactory::FlexibleInterpVar* pw = dynamic_cast<const RooStats::HistFactory::FlexibleInterpVar*>(pw_);
if (!pw) {
std::cout<<"Skipping object '"<<pw_->GetName()<<"': not a FlexibleInterpVar"<<std::endl;
continue;
}
processFlexibleInterpVar(pw, pattern, vetoPattern, scale);
}
}
return true;
}
};
namespace {
bool available = TSStatisticsManager::registerAction(new ScaleSystematics(),"ScaleSystematics");
}
}