Logo ROOT  
Reference Guide
StandardHistFactoryPlotsWithCategories.C File Reference

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`
float * q
Definition: THbookFile.cxx:87
␛[1mRooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby␛[0m
Copyright (C) 2000-2013 NIKHEF, University of California & Stanford University
All rights reserved, please read http://roofit.sourceforge.net/license.txt
1) 0x55b213635720 RooRealVar:: alpha_syst2 = 0 L(-5 - 5) "alpha_syst2"
2) 0x55b21363b840 RooRealVar:: alpha_syst3 = 0 L(-5 - 5) "alpha_syst3"
3) 0x55b211d5c900 RooRealVar:: gamma_stat_channel1_bin_0 = 1 +/- 0.05 L(0 - 1.25) "gamma_stat_channel1_bin_0"
4) 0x55b211db2b80 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:InputArguments -- The formula select claims to use the variables (channelCat,obs_x_channel1) but only (channelCat) seem to be in use.
inputs: channelCat==channelCat::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)
[#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:InputArguments -- The formula select claims to use the variables (channelCat,obs_x_channel1) but only (channelCat) seem to be in use.
inputs: channelCat==channelCat::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:InputArguments -- The formula select claims to use the variables (channelCat,obs_x_channel1) but only (channelCat) seem to be in use.
inputs: channelCat==channelCat::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:InputArguments -- The formula select claims to use the variables (channelCat,obs_x_channel1) but only (channelCat) seem to be in use.
inputs: channelCat==channelCat::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)
#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 "RooCategory.h"
using namespace RooFit;
using namespace RooStats;
using namespace std;
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) {
#ifdef _WIN32
cout << "HistFactory file cannot be generated on Windows - exit" << endl;
return;
#endif
// 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 byt 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();
TList *list = new TList();
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());
TIterator *iter = channelCat->typeIterator();
RooCatType *tt = NULL;
tt = (RooCatType *)iter->Next();
RooAbsPdf *pdftmp = ((RooSimultaneous *)mc->GetPdf())->getPdf(tt->GetName());
RooArgSet *obstmp = pdftmp->getObservables(*mc->GetObservables());
obs = ((RooRealVar *)obstmp->first());
RooPlot *frame = obs->frame();
cout << Form("%s==%s::%s", channelCat->GetName(), channelCat->GetName(), tt->GetName()) << endl;
cout << tt->GetName() << " " << channelCat->getLabel() << endl;
data->plotOn(frame, MarkerSize(1),
Cut(Form("%s==%s::%s", channelCat->GetName(), channelCat->GetName(), tt->GetName())),
Double_t normCount =
data->sumEntries(Form("%s==%s::%s", channelCat->GetName(), channelCat->GetName(), tt->GetName()));
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) {
TIterator *it = mc->GetNuisanceParameters()->createIterator();
RooRealVar *var = NULL;
while ((var = (RooRealVar *)it->Next()) != NULL) {
RooPlot *frame = obs->frame();
frame->SetYTitle(var->GetName());
data->plotOn(frame, MarkerSize(1));
var->setVal(0);
mc->GetPdf()->plotOn(frame, LineWidth(1.));
var->setVal(1);
mc->GetPdf()->plotOn(frame, LineColor(kRed), LineStyle(kDashed), LineWidth(1));
var->setVal(-1);
mc->GetPdf()->plotOn(frame, LineColor(kGreen), LineStyle(kDashed), LineWidth(1));
list->Add(frame);
var->setVal(0);
}
} else {
RooCategory *channelCat = (RooCategory *)(&simPdf->indexCat());
// TIterator* iter = simPdf->indexCat().typeIterator() ;
TIterator *iter = channelCat->typeIterator();
RooCatType *tt = NULL;
while (nPlots < nPlotsMax && (tt = (RooCatType *)iter->Next())) {
cout << "on type " << tt->GetName() << " " << endl;
// Get pdf associated with state from simpdf
RooAbsPdf *pdftmp = simPdf->getPdf(tt->GetName());
// Generate observables defined by the pdf associated with this state
RooArgSet *obstmp = pdftmp->getObservables(*mc->GetObservables());
// obstmp->Print();
obs = ((RooRealVar *)obstmp->first());
TIterator *it = mc->GetNuisanceParameters()->createIterator();
RooRealVar *var = NULL;
while (nPlots < nPlotsMax && (var = (RooRealVar *)it->Next())) {
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(), tt->GetName()) << endl;
cout << tt->GetName() << " " << channelCat->getLabel() << endl;
data->plotOn(frame, MarkerSize(1),
Cut(Form("%s==%s::%s", channelCat->GetName(), channelCat->GetName(), tt->GetName())),
Double_t normCount =
data->sumEntries(Form("%s==%s::%s", channelCat->GetName(), channelCat->GetName(), tt->GetName()));
if (strcmp(var->GetName(), "Lumi") == 0) {
cout << "working on lumi" << endl;
var->setVal(w->var("nominalLumi")->getVal());
var->Print();
} else {
var->setVal(0);
}
// w->allVars().Print("v");
// mc->GetNuisanceParameters()->Print("v");
// pdftmp->plotOn(frame,LineWidth(2.));
// mc->GetPdf()->plotOn(frame,LineWidth(2.),Slice(*channelCat,tt->GetName()),ProjWData(*data));
// pdftmp->plotOn(frame,LineWidth(2.),Slice(*channelCat,tt->GetName()),ProjWData(*data));
normCount = pdftmp->expectedEvents(*obs);
pdftmp->plotOn(frame, LineWidth(2.), Normalization(normCount, RooAbsReal::NumEvent));
if (strcmp(var->GetName(), "Lumi") == 0) {
cout << "working on lumi" << endl;
var->setVal(w->var("nominalLumi")->getVal() + 0.05);
var->Print();
} else {
var->setVal(nSigmaToVary);
}
// pdftmp->plotOn(frame,LineColor(kRed),LineStyle(kDashed),LineWidth(2));
// mc->GetPdf()->plotOn(frame,LineColor(kRed),LineStyle(kDashed),LineWidth(2.),Slice(*channelCat,tt->GetName()),ProjWData(*data));
// pdftmp->plotOn(frame,LineColor(kRed),LineStyle(kDashed),LineWidth(2.),Slice(*channelCat,tt->GetName()),ProjWData(*data));
normCount = pdftmp->expectedEvents(*obs);
pdftmp->plotOn(frame, LineWidth(2.), LineColor(kRed), LineStyle(kDashed),
if (strcmp(var->GetName(), "Lumi") == 0) {
cout << "working on lumi" << endl;
var->setVal(w->var("nominalLumi")->getVal() - 0.05);
var->Print();
} else {
var->setVal(-nSigmaToVary);
}
// pdftmp->plotOn(frame,LineColor(kGreen),LineStyle(kDashed),LineWidth(2));
// mc->GetPdf()->plotOn(frame,LineColor(kGreen),LineStyle(kDashed),LineWidth(2),Slice(*channelCat,tt->GetName()),ProjWData(*data));
// pdftmp->plotOn(frame,LineColor(kGreen),LineStyle(kDashed),LineWidth(2),Slice(*channelCat,tt->GetName()),ProjWData(*data));
normCount = pdftmp->expectedEvents(*obs);
// set them back to normal
if (strcmp(var->GetName(), "Lumi") == 0) {
cout << "working on lumi" << endl;
var->setVal(w->var("nominalLumi")->getVal());
var->Print();
} else {
var->setVal(0);
}
list->Add(frame);
// quit making plots
++nPlots;
frame->Draw();
c2->SaveAs(Form("%s_%s_%s.pdf", tt->GetName(), obs->GetName(), var->GetName()));
delete c2;
}
}
}
// -------------------------------------------------------
// now make plots
TCanvas *c1 = new TCanvas("c1", "ProfileInspectorDemo", 800, 200);
if (list->GetSize() > 4) {
double n = list->GetSize();
int nx = (int)sqrt(n);
int ny = TMath::CeilNint(n / nx);
c1->Divide(ny, nx);
} else
c1->Divide(list->GetSize());
for (int i = 0; i < list->GetSize(); ++i) {
c1->cd(i + 1);
list->At(i)->Draw();
}
}
double Double_t
Definition: RtypesCore.h:57
@ kRed
Definition: Rtypes.h:64
@ kGreen
Definition: Rtypes.h:64
@ kDashed
Definition: TAttLine.h:48
double sqrt(double)
#define gROOT
Definition: TROOT.h:406
char * Form(const char *fmt,...)
R__EXTERN TSystem * gSystem
Definition: TSystem.h:556
RooArgSet * getObservables(const RooArgSet &set, Bool_t valueOnly=kTRUE) const
Return the observables of this pdf given a set of observables.
Definition: RooAbsArg.h:276
virtual void Print(Option_t *options=0) const
Print the object to the defaultPrintStream().
Definition: RooAbsArg.h:302
TIterator * typeIterator() const
const char * getLabel() const
Retrieve current label. Use getCurrentLabel() for more clarity.
RooAbsArg * first() const
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:44
virtual const RooArgSet * get() const
Definition: RooAbsData.h:87
virtual Double_t sumEntries() const =0
virtual RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none()) const
Calls RooPlot* plotOn(RooPlot* frame, const RooLinkedList& cmdList) const ;.
Definition: RooAbsData.cxx:546
virtual Double_t expectedEvents(const RooArgSet *nset) const
Return expected number of events from this p.d.f for use in extended likelihood calculations.
Definition: RooAbsPdf.cxx:3283
virtual RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none(), const RooCmdArg &arg9=RooCmdArg::none(), const RooCmdArg &arg10=RooCmdArg::none()) const
Helper calling plotOn(RooPlot*, RooLinkedList&) const.
Definition: RooAbsPdf.h:118
RooPlot * frame(const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none()) const
Create a new RooPlot on the heap with a drawing frame initialized for this object,...
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:90
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
RooCatType is an auxilary class for RooAbsCategory and defines a a single category state.
RooCategory is an object to represent discrete states.
Definition: RooCategory.h:23
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
A RooPlot is a plot frame and a container for graphics objects within that frame.
Definition: RooPlot.h:44
void SetName(const char *name)
Set the name of the RooPlot to 'name'.
Definition: RooPlot.cxx:1236
void SetYTitle(const char *title)
Definition: RooPlot.cxx:1314
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:712
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:261
RooSimultaneous facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset.
const RooAbsCategoryLValue & indexCat() const
RooAbsPdf * getPdf(const char *catName) const
Return the p.d.f associated with the given index category name.
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
RooAbsData * data(const char *name) const
Retrieve dataset (binned or unbinned) with given name. A null pointer is returned if not found.
void Print(Option_t *opts=0) const
Print contents of the workspace.
RooRealVar * var(const char *name) const
Retrieve real-valued variable (RooRealVar) with given name. A null pointer is returned if not found.
TObject * obj(const char *name) const
Return any type of object (RooAbsArg, RooAbsData or generic object) with given name)
The Canvas class.
Definition: TCanvas.h:27
virtual Int_t GetSize() const
Return the capacity of the collection, i.e.
Definition: TCollection.h:182
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:53
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:3942
Iterator abstract base class.
Definition: TIterator.h:30
virtual TObject * Next()=0
A doubly linked list.
Definition: TList.h:44
virtual void Add(TObject *obj)
Definition: TList.h:87
virtual TObject * At(Int_t idx) const
Returns the object at position idx. Returns 0 if idx is out of range.
Definition: TList.cxx:356
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition: TObject.cxx:195
virtual Bool_t AccessPathName(const char *path, EAccessMode mode=kFileExists)
Returns FALSE if one can access a file using the specified access mode.
Definition: TSystem.cxx:1291
RooCmdArg DataError(Int_t)
RooCmdArg LineWidth(Width_t width)
RooCmdArg MarkerSize(Size_t size)
RooCmdArg Cut(const char *cutSpec)
RooCmdArg LineColor(Color_t color)
RooCmdArg Normalization(Double_t scaleFactor)
RooCmdArg LineStyle(Style_t style)
return c1
Definition: legend1.C:41
const Int_t n
Definition: legend1.C:16
return c2
Definition: legend2.C:14
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Namespace for the RooStats classes.
Definition: Asimov.h:19
Int_t CeilNint(Double_t x)
Definition: TMath.h:689
Definition: file.py:1
auto * tt
Definition: textangle.C:16
Author
Kyle Cranmer

Definition in file StandardHistFactoryPlotsWithCategories.C.