Logo ROOT  

Reference Guide
Loading...
Searching...
No Matches
StandardHistFactoryPlotsWithCategories.C File Reference
Tutorials » RooFit packages tutorials » RooStats

Detailed Description

View in nbviewer Open in SWAN
StandardHistFactoryPlotsWithCategories

This is a standard demo that can be used with any ROOT file prepared in the standard way. You specify:

  • name for input ROOT file
  • name of workspace inside ROOT file that holds model and data
  • name of ModelConfig that specifies details for calculator tools
  • name of dataset

With default parameters the macro will attempt to run the standard hist2workspace example and read the ROOT file that it produces.

The macro will scan through all the categories in a simPdf find the corresponding observable. For each category, it will loop through each of the nuisance parameters and plot

  • the data
  • the nominal model (blue)
  • the +Nsigma (red)
  • the -Nsigma (green)

You can specify how many sigma to vary by changing nSigmaToVary. You can also change the signal rate by changing muVal.

The script produces a lot plots, you can merge them by doing:

gs -q -dNOPAUSE -dBATCH -sDEVICE=pdfwrite -sOutputFile=merged.pdf `ls *pdf`
q
float * q
Definition THbookFile.cxx:89
1) 0x559558e8b780 RooRealVar:: alpha_syst2 = 0 +/- 1 L(-5 - 5) "alpha_syst2"
2) 0x559559d19340 RooRealVar:: alpha_syst3 = 0 +/- 1 L(-5 - 5) "alpha_syst3"
3) 0x559559ca86d0 RooRealVar:: gamma_stat_channel1_bin_0 = 1 +/- 0.05 L(0 - 1.25) "gamma_stat_channel1_bin_0"
4) 0x559558c13610 RooRealVar:: gamma_stat_channel1_bin_1 = 1 +/- 0.1 L(0 - 1.5) "gamma_stat_channel1_bin_1"
check expectedData by category
Is a simultaneous PDF
on type channel1
channelCat==channelCat::channel1
channel1 channel1
[#1] INFO:Plotting -- RooTreeData::plotOn: plotting 234 events out of 234 total events
[#1] INFO:NumericIntegration -- RooRealIntegral::init(channel1_model_Int[obs_x_channel1]) using numeric integrator RooBinIntegrator to calculate Int(obs_x_channel1)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(channel1_model_Int[obs_x_channel1]) using numeric integrator RooBinIntegrator to calculate Int(obs_x_channel1)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(channel1_model_Int[obs_x_channel1]) using numeric integrator RooBinIntegrator to calculate Int(obs_x_channel1)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(channel1_model_Int[obs_x_channel1]) using numeric integrator RooBinIntegrator to calculate Int(obs_x_channel1)
channelCat==channelCat::channel1
channel1 channel1
[#1] INFO:Plotting -- RooTreeData::plotOn: plotting 234 events out of 234 total events
[#1] INFO:NumericIntegration -- RooRealIntegral::init(channel1_model_Int[obs_x_channel1]) using numeric integrator RooBinIntegrator to calculate Int(obs_x_channel1)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(channel1_model_Int[obs_x_channel1]) using numeric integrator RooBinIntegrator to calculate Int(obs_x_channel1)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(channel1_model_Int[obs_x_channel1]) using numeric integrator RooBinIntegrator to calculate Int(obs_x_channel1)
channelCat==channelCat::channel1
channel1 channel1
[#1] INFO:Plotting -- RooTreeData::plotOn: plotting 234 events out of 234 total events
[#1] INFO:NumericIntegration -- RooRealIntegral::init(channel1_model_Int[obs_x_channel1]) using numeric integrator RooBinIntegrator to calculate Int(obs_x_channel1)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(channel1_model_Int[obs_x_channel1]) using numeric integrator RooBinIntegrator to calculate Int(obs_x_channel1)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(channel1_model_Int[obs_x_channel1]) using numeric integrator RooBinIntegrator to calculate Int(obs_x_channel1)
channelCat==channelCat::channel1
channel1 channel1
[#1] INFO:Plotting -- RooTreeData::plotOn: plotting 234 events out of 234 total events
[#1] INFO:NumericIntegration -- RooRealIntegral::init(channel1_model_Int[obs_x_channel1]) using numeric integrator RooBinIntegrator to calculate Int(obs_x_channel1)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(channel1_model_Int[obs_x_channel1]) using numeric integrator RooBinIntegrator to calculate Int(obs_x_channel1)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(channel1_model_Int[obs_x_channel1]) using numeric integrator RooBinIntegrator to calculate Int(obs_x_channel1)
#include "TFile.h"
#include "TROOT.h"
#include "TCanvas.h"
#include "TList.h"
#include "TMath.h"
#include "TSystem.h"
#include "RooWorkspace.h"
#include "RooAbsData.h"
#include "RooRealVar.h"
#include "RooPlot.h"
#include "RooSimultaneous.h"
#include "RooCategory.h"
#include "RooStats/ModelConfig.h"
#include "RooStats/ProfileInspector.h"
using namespace RooFit;
using namespace RooStats;
using std::cout, std::endl;
void StandardHistFactoryPlotsWithCategories(const char *infile = "", const char *workspaceName = "combined",
const char *modelConfigName = "ModelConfig",
const char *dataName = "obsData")
{
double nSigmaToVary = 5.;
double muVal = 0;
bool doFit = false;
// -------------------------------------------------------
// First part is just to access a user-defined file
// or create the standard example file if it doesn't exist
const char *filename = "";
if (!strcmp(infile, "")) {
filename = "results/example_combined_GaussExample_model.root";
bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code
// if file does not exists generate with histfactory
if (!fileExist) {
// Normally this would be run on the command line
cout << "will run standard hist2workspace example" << endl;
gROOT->ProcessLine(".! prepareHistFactory .");
gROOT->ProcessLine(".! hist2workspace config/example.xml");
cout << "\n\n---------------------" << endl;
cout << "Done creating example input" << endl;
cout << "---------------------\n\n" << endl;
}
} else
filename = infile;
// Try to open the file
TFile *file = TFile::Open(filename);
// if input file was specified but not found, quit
if (!file) {
cout << "StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
return;
}
// -------------------------------------------------------
// Tutorial starts here
// -------------------------------------------------------
// get the workspace out of the file
RooWorkspace *w = (RooWorkspace *)file->Get(workspaceName);
if (!w) {
cout << "workspace not found" << endl;
return;
}
// get the modelConfig out of the file
ModelConfig *mc = (ModelConfig *)w->obj(modelConfigName);
// get the modelConfig out of the file
RooAbsData *data = w->data(dataName);
// make sure ingredients are found
if (!data || !mc) {
w->Print();
cout << "data or ModelConfig was not found" << endl;
return;
}
// -------------------------------------------------------
// now use the profile inspector
RooRealVar *obs = (RooRealVar *)mc->GetObservables()->first();
std::vector<RooPlot *> frameList;
RooRealVar *firstPOI = dynamic_cast<RooRealVar *>(mc->GetParametersOfInterest()->first());
firstPOI->setVal(muVal);
// firstPOI->setConstant();
if (doFit) {
mc->GetPdf()->fitTo(*data);
}
// -------------------------------------------------------
mc->GetNuisanceParameters()->Print("v");
int nPlotsMax = 1000;
cout << " check expectedData by category" << endl;
RooDataSet *simData = NULL;
RooSimultaneous *simPdf = NULL;
if (strcmp(mc->GetPdf()->ClassName(), "RooSimultaneous") == 0) {
cout << "Is a simultaneous PDF" << endl;
simPdf = (RooSimultaneous *)(mc->GetPdf());
} else {
cout << "Is not a simultaneous PDF" << endl;
}
if (doFit) {
RooCategory *channelCat = (RooCategory *)(&simPdf->indexCat());
auto const& catName = channelCat->begin()->first;
RooAbsPdf *pdftmp = ((RooSimultaneous *)mc->GetPdf())->getPdf(catName.c_str());
std::unique_ptr<RooArgSet> obstmp{pdftmp->getObservables(*mc->GetObservables())};
obs = ((RooRealVar *)obstmp->first());
RooPlot *frame = obs->frame();
cout << Form("%s==%s::%s", channelCat->GetName(), channelCat->GetName(), catName.c_str()) << endl;
cout << catName << " " << channelCat->getLabel() << endl;
data->plotOn(frame, MarkerSize(1),
Cut(Form("%s==%s::%s", channelCat->GetName(), channelCat->GetName(), catName.c_str())),
DataError(RooAbsData::None));
Double_t normCount =
data->sumEntries(Form("%s==%s::%s", channelCat->GetName(), channelCat->GetName(), catName.c_str()));
pdftmp->plotOn(frame, LineWidth(2.), Normalization(normCount, RooAbsReal::NumEvent));
frame->Draw();
cout << "expected events = " << mc->GetPdf()->expectedEvents(*data->get()) << endl;
return;
}
int nPlots = 0;
if (!simPdf) {
for (auto *var : static_range_cast<RooRealVar *>(*mc->GetNuisanceParameters())) {
RooPlot *frame = obs->frame();
frame->SetYTitle(var->GetName());
data->plotOn(frame, MarkerSize(1));
const double value = var->getVal();
mc->GetPdf()->plotOn(frame, LineWidth(1.));
var->setVal(value + var->getError());
mc->GetPdf()->plotOn(frame, LineColor(kRed), LineStyle(kDashed), LineWidth(1));
var->setVal(value - var->getError());
mc->GetPdf()->plotOn(frame, LineColor(kGreen), LineStyle(kDashed), LineWidth(1));
frameList.push_back(frame);
var->setVal(value);
}
} else {
RooCategory *channelCat = (RooCategory *)(&simPdf->indexCat());
for (auto const& tt : *channelCat) {
if (nPlots == nPlotsMax) {
break;
}
auto const& catName = tt.first;
cout << "on type " << catName << " " << endl;
// Get pdf associated with state from simpdf
RooAbsPdf *pdftmp = simPdf->getPdf(catName.c_str());
// Generate observables defined by the pdf associated with this state
std::unique_ptr<RooArgSet> obstmp{pdftmp->getObservables(*mc->GetObservables())};
// obstmp->Print();
obs = ((RooRealVar *)obstmp->first());
for (auto *var : static_range_cast<RooRealVar*>(*mc->GetNuisanceParameters())) {
if (nPlots == nPlotsMax) break;
TCanvas *c2 = new TCanvas("c2");
RooPlot *frame = obs->frame();
frame->SetName(Form("frame%d", nPlots));
frame->SetYTitle(var->GetName());
cout << Form("%s==%s::%s", channelCat->GetName(), channelCat->GetName(), catName.c_str()) << endl;
cout << catName << " " << channelCat->getLabel() << endl;
data->plotOn(frame, MarkerSize(1),
Cut(Form("%s==%s::%s", channelCat->GetName(), channelCat->GetName(), catName.c_str())),
DataError(RooAbsData::None));
Double_t normCount =
data->sumEntries(Form("%s==%s::%s", channelCat->GetName(), channelCat->GetName(), catName.c_str()));
// remember the nominal value
const double value = var->getVal();
// w->allVars().Print("v");
// mc->GetNuisanceParameters()->Print("v");
// pdftmp->plotOn(frame,LineWidth(2.));
// mc->GetPdf()->plotOn(frame,LineWidth(2.),Slice(*channelCat,catName.c_str()),ProjWData(*data));
// pdftmp->plotOn(frame,LineWidth(2.),Slice(*channelCat,catName.c_str()),ProjWData(*data));
normCount = pdftmp->expectedEvents(*obs);
pdftmp->plotOn(frame, LineWidth(2.), Normalization(normCount, RooAbsReal::NumEvent)); // nominal
var->setVal(value + nSigmaToVary * var->getError());
// pdftmp->plotOn(frame,LineColor(kRed),LineStyle(kDashed),LineWidth(2));
// mc->GetPdf()->plotOn(frame,LineColor(kRed),LineStyle(kDashed),LineWidth(2.),Slice(*channelCat,catName.c_str()),ProjWData(*data));
// pdftmp->plotOn(frame,LineColor(kRed),LineStyle(kDashed),LineWidth(2.),Slice(*channelCat,catName.c_str()),ProjWData(*data));
normCount = pdftmp->expectedEvents(*obs);
pdftmp->plotOn(frame, LineWidth(2.), LineColor(kRed), LineStyle(kDashed),
Normalization(normCount, RooAbsReal::NumEvent)); // +n sigma
var->setVal(value - nSigmaToVary * var->getError());
// pdftmp->plotOn(frame,LineColor(kGreen),LineStyle(kDashed),LineWidth(2));
// mc->GetPdf()->plotOn(frame,LineColor(kGreen),LineStyle(kDashed),LineWidth(2),Slice(*channelCat,catName.c_str()),ProjWData(*data));
// pdftmp->plotOn(frame,LineColor(kGreen),LineStyle(kDashed),LineWidth(2),Slice(*channelCat,catName.c_str()),ProjWData(*data));
normCount = pdftmp->expectedEvents(*obs);
pdftmp->plotOn(frame, LineWidth(2.), LineColor(kGreen), LineStyle(kDashed),
Normalization(normCount, RooAbsReal::NumEvent)); // -n sigma
// set them back to normal
var->setVal(value);
frameList.push_back(frame);
// quit making plots
++nPlots;
frame->Draw();
c2->SaveAs(Form("%s_%s_%s.pdf", catName.c_str(), obs->GetName(), var->GetName()));
delete c2;
}
}
}
// -------------------------------------------------------
// now make plots
TCanvas *c1 = new TCanvas("c1", "ProfileInspectorDemo", 800, 200);
int nFrames = frameList.size();
if (nFrames > 4) {
int nx = (int)sqrt(nFrames);
int ny = TMath::CeilNint(nFrames / nx);
nx = TMath::CeilNint(sqrt(nFrames));
c1->Divide(ny, nx);
} else
c1->Divide(nFrames);
for (int i = 0; i < nFrames; ++i) {
c1->cd(i + 1);
frameList[i]->Draw();
}
}
static_range_cast
ROOT::RRangeCast< T, false, Range_t > static_range_cast(Range_t &&coll)
Definition RooAbsCollection.h:45
Double_t
double Double_t
Double 8 bytes.
Definition RtypesCore.h:73
kRed
@ kRed
Definition Rtypes.h:67
kGreen
@ kGreen
Definition Rtypes.h:67
kDashed
@ kDashed
Definition TAttLine.h:54
gROOT
#define gROOT
Definition TROOT.h:417
Form
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition TString.cxx:2496
gSystem
externTSystem * gSystem
Definition TSystem.h:582
RooAbsArg::getObservables
RooFit::OwningPtr< RooArgSet > getObservables(const RooArgSet &set, bool valueOnly=true) const
Given a set of possible observables, return the observables that this PDF depends on.
Definition RooAbsArg.cxx:596
RooAbsCategory::getLabel
const char * getLabel() const
Retrieve current label. Use getCurrentLabel() for more clarity.
Definition RooAbsCategory.h:136
RooAbsCategory::begin
std::map< std::string, value_type >::const_iterator begin() const
Iterator for category state names. Points to pairs of index and name.
Definition RooAbsCategory.h:92
RooAbsCollection::first
RooAbsArg * first() const
Definition RooAbsCollection.h:296
RooAbsCollection::Print
void Print(Option_t *options=nullptr) const override
This method must be overridden when a class wants to print itself.
Definition RooAbsCollection.h:320
RooAbsData
Abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:56
RooAbsPdf
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:32
RooAbsPdf::plotOn
RooPlot * plotOn(RooPlot *frame, 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={}, const RooCmdArg &arg10={}) const override
Helper calling plotOn(RooPlot*, RooLinkedList&) const.
Definition RooAbsPdf.h:116
RooAbsPdf::expectedEvents
virtual double expectedEvents(const RooArgSet *nset) const
Return expected number of events to be used in calculation of extended likelihood.
Definition RooAbsPdf.cxx:2398
RooAbsPdf::fitTo
RooFit::OwningPtr< RooFitResult > fitTo(RooAbsData &data, CmdArgs_t const &... cmdArgs)
Fit PDF to given dataset.
Definition RooAbsPdf.h:149
RooAbsRealLValue::frame
RooPlot * frame(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
Create a new RooPlot on the heap with a drawing frame initialized for this object,...
Definition RooAbsRealLValue.cxx:199
RooCategory
Object to represent discrete states.
Definition RooCategory.h:28
RooDataSet
Container class to hold unbinned data.
Definition RooDataSet.h:32
RooPlot
Plot frame and a container for graphics objects within that frame.
Definition RooPlot.h:43
RooPlot::SetName
void SetName(const char *name) override
Set the name of the RooPlot to 'name'.
Definition RooPlot.cxx:1187
RooPlot::Draw
void Draw(Option_t *options=nullptr) override
Draw this plot and all of the elements it contains.
Definition RooPlot.cxx:597
RooPlot::SetYTitle
void SetYTitle(const char *title)
Definition RooPlot.cxx:1302
RooRealVar
Variable that can be changed from the outside.
Definition RooRealVar.h:37
RooRealVar::setVal
void setVal(double value) override
Set value of variable to 'value'.
Definition RooRealVar.cxx:284
RooSimultaneous
Facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset.
Definition RooSimultaneous.h:41
RooSimultaneous::getPdf
RooAbsPdf * getPdf(RooStringView catName) const
Return the p.d.f associated with the given index category name.
Definition RooSimultaneous.cxx:376
RooSimultaneous::indexCat
const RooAbsCategoryLValue & indexCat() const
Definition RooSimultaneous.h:91
RooStats::ModelConfig::GetParametersOfInterest
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return nullptr if not existing)
Definition ModelConfig.h:273
RooStats::ModelConfig::GetNuisanceParameters
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return nullptr if not existing)
Definition ModelConfig.h:276
RooStats::ModelConfig::GetObservables
const RooArgSet * GetObservables() const
get RooArgSet for observables (return nullptr if not existing)
Definition ModelConfig.h:285
RooStats::ModelConfig::GetPdf
RooAbsPdf * GetPdf() const
get model PDF (return nullptr if pdf has not been specified or does not exist)
Definition ModelConfig.h:270
RooWorkspace
Persistable container for RooFit projects.
Definition RooWorkspace.h:43
TCanvas
The Canvas class.
Definition TCanvas.h:23
TDirectoryFile::Get
TObject * Get(const char *namecycle) override
Return pointer to object identified by namecycle.
TFile
A file, usually with extension .root, that stores data and code in the form of serialized objects in ...
Definition TFile.h:130
TFile::Open
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=ROOT::RCompressionSetting::EDefaults::kUseCompiledDefault, Int_t netopt=0)
Create / open a file.
Definition TFile.cxx:3787
TNamed::GetName
const char * GetName() const override
Returns name of object.
Definition TNamed.h:49
TObject::ClassName
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition TObject.cxx:227
RooFit::DataError
RooCmdArg DataError(Int_t)
Definition RooGlobalFunc.cxx:455
RooFit::LineColor
RooCmdArg LineColor(TColorNumber color)
Definition RooGlobalFunc.cxx:229
RooFit::LineWidth
RooCmdArg LineWidth(Width_t width)
Definition RooGlobalFunc.cxx:241
RooFit::Normalization
RooCmdArg Normalization(double scaleFactor)
Definition RooGlobalFunc.cxx:209
RooFit::MarkerSize
RooCmdArg MarkerSize(Size_t size)
Definition RooGlobalFunc.cxx:337
RooFit::Cut
RooCmdArg Cut(const char *cutSpec)
Definition RooGlobalFunc.cxx:309
RooFit::LineStyle
RooCmdArg LineStyle(Style_t style)
Definition RooGlobalFunc.cxx:233
c1
return c1
Definition legend1.C:41
c2
return c2
Definition legend2.C:14
ROOT::Math::sqrt
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
Definition UnaryOperators.h:281
RooFit
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition CodegenImpl.h:72
RooFit::ModelConfig
RooStats::ModelConfig ModelConfig
Definition ModelConfig.h:377
RooStats
Namespace for the RooStats classes.
Definition CodegenImpl.h:66
TMath::CeilNint
Int_t CeilNint(Double_t x)
Returns the nearest integer of TMath::Ceil(x).
Definition TMath.h:685
rf503_wspaceread.data
data
Definition rf503_wspaceread.py:37
rf503_wspaceread.w
w
Definition rf503_wspaceread.py:29
tt
auto * tt
Definition textangle.C:16
Author
Kyle Cranmer

Definition in file StandardHistFactoryPlotsWithCategories.C.