Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
StandardHistFactoryPlotsWithCategories.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_roostats
3/// \notebook -js
4/// StandardHistFactoryPlotsWithCategories
5///
6/// This is a standard demo that can be used with any ROOT file
7/// prepared in the standard way. You specify:
8/// - name for input ROOT file
9/// - name of workspace inside ROOT file that holds model and data
10/// - name of ModelConfig that specifies details for calculator tools
11/// - name of dataset
12///
13/// With default parameters the macro will attempt to run the
14/// standard hist2workspace example and read the ROOT file
15/// that it produces.
16///
17/// The macro will scan through all the categories in a simPdf find the corresponding
18/// observable. For each category, it will loop through each of the nuisance parameters
19/// and plot
20/// - the data
21/// - the nominal model (blue)
22/// - the +Nsigma (red)
23/// - the -Nsigma (green)
24///
25/// You can specify how many sigma to vary by changing nSigmaToVary.
26/// You can also change the signal rate by changing muVal.
27///
28/// The script produces a lot plots, you can merge them by doing:
29/// ~~~{.cpp}
30/// gs -q -dNOPAUSE -dBATCH -sDEVICE=pdfwrite -sOutputFile=merged.pdf `ls *pdf`
31/// ~~~
32///
33/// \macro_image
34/// \macro_output
35/// \macro_code
36///
37/// \author Kyle Cranmer
38
39#include "TFile.h"
40#include "TROOT.h"
41#include "TCanvas.h"
42#include "TList.h"
43#include "TMath.h"
44#include "TSystem.h"
45#include "RooWorkspace.h"
46#include "RooAbsData.h"
47#include "RooRealVar.h"
48#include "RooPlot.h"
49#include "RooSimultaneous.h"
50#include "RooCategory.h"
51
54
55using namespace RooFit;
56using namespace RooStats;
57using namespace std;
58
59void StandardHistFactoryPlotsWithCategories(const char *infile = "", const char *workspaceName = "combined",
60 const char *modelConfigName = "ModelConfig",
61 const char *dataName = "obsData")
62{
63
64 double nSigmaToVary = 5.;
65 double muVal = 0;
66 bool doFit = false;
67
68 // -------------------------------------------------------
69 // First part is just to access a user-defined file
70 // or create the standard example file if it doesn't exist
71 const char *filename = "";
72 if (!strcmp(infile, "")) {
73 filename = "results/example_combined_GaussExample_model.root";
74 bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code
75 // if file does not exists generate with histfactory
76 if (!fileExist) {
77#ifdef _WIN32
78 cout << "HistFactory file cannot be generated on Windows - exit" << endl;
79 return;
80#endif
81 // Normally this would be run on the command line
82 cout << "will run standard hist2workspace example" << endl;
83 gROOT->ProcessLine(".! prepareHistFactory .");
84 gROOT->ProcessLine(".! hist2workspace config/example.xml");
85 cout << "\n\n---------------------" << endl;
86 cout << "Done creating example input" << endl;
87 cout << "---------------------\n\n" << endl;
88 }
89
90 } else
91 filename = infile;
92
93 // Try to open the file
94 TFile *file = TFile::Open(filename);
95
96 // if input file was specified byt not found, quit
97 if (!file) {
98 cout << "StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
99 return;
100 }
101
102 // -------------------------------------------------------
103 // Tutorial starts here
104 // -------------------------------------------------------
105
106 // get the workspace out of the file
107 RooWorkspace *w = (RooWorkspace *)file->Get(workspaceName);
108 if (!w) {
109 cout << "workspace not found" << endl;
110 return;
111 }
112
113 // get the modelConfig out of the file
114 ModelConfig *mc = (ModelConfig *)w->obj(modelConfigName);
115
116 // get the modelConfig out of the file
117 RooAbsData *data = w->data(dataName);
118
119 // make sure ingredients are found
120 if (!data || !mc) {
121 w->Print();
122 cout << "data or ModelConfig was not found" << endl;
123 return;
124 }
125
126 // -------------------------------------------------------
127 // now use the profile inspector
128
129 RooRealVar *obs = (RooRealVar *)mc->GetObservables()->first();
130 TList *list = new TList();
131
132 RooRealVar *firstPOI = dynamic_cast<RooRealVar *>(mc->GetParametersOfInterest()->first());
133
134 firstPOI->setVal(muVal);
135 // firstPOI->setConstant();
136 if (doFit) {
137 mc->GetPdf()->fitTo(*data);
138 }
139
140 // -------------------------------------------------------
141
142 mc->GetNuisanceParameters()->Print("v");
143 int nPlotsMax = 1000;
144 cout << " check expectedData by category" << endl;
145 RooDataSet *simData = NULL;
146 RooSimultaneous *simPdf = NULL;
147 if (strcmp(mc->GetPdf()->ClassName(), "RooSimultaneous") == 0) {
148 cout << "Is a simultaneous PDF" << endl;
149 simPdf = (RooSimultaneous *)(mc->GetPdf());
150 } else {
151 cout << "Is not a simultaneous PDF" << endl;
152 }
153
154 if (doFit) {
155 RooCategory *channelCat = (RooCategory *)(&simPdf->indexCat());
156 TIterator *iter = channelCat->typeIterator();
157 RooCatType *tt = NULL;
158 tt = (RooCatType *)iter->Next();
159 RooAbsPdf *pdftmp = ((RooSimultaneous *)mc->GetPdf())->getPdf(tt->GetName());
160 RooArgSet *obstmp = pdftmp->getObservables(*mc->GetObservables());
161 obs = ((RooRealVar *)obstmp->first());
162 RooPlot *frame = obs->frame();
163 cout << Form("%s==%s::%s", channelCat->GetName(), channelCat->GetName(), tt->GetName()) << endl;
164 cout << tt->GetName() << " " << channelCat->getLabel() << endl;
165 data->plotOn(frame, MarkerSize(1),
166 Cut(Form("%s==%s::%s", channelCat->GetName(), channelCat->GetName(), tt->GetName())),
168
169 Double_t normCount =
170 data->sumEntries(Form("%s==%s::%s", channelCat->GetName(), channelCat->GetName(), tt->GetName()));
171
172 pdftmp->plotOn(frame, LineWidth(2.), Normalization(normCount, RooAbsReal::NumEvent));
173 frame->Draw();
174 cout << "expected events = " << mc->GetPdf()->expectedEvents(*data->get()) << endl;
175 return;
176 }
177
178 int nPlots = 0;
179 if (!simPdf) {
180
182 RooRealVar *var = NULL;
183 while ((var = (RooRealVar *)it.Next()) != NULL) {
184 RooPlot *frame = obs->frame();
185 frame->SetYTitle(var->GetName());
186 data->plotOn(frame, MarkerSize(1));
187 var->setVal(0);
188 mc->GetPdf()->plotOn(frame, LineWidth(1.));
189 var->setVal(1);
190 mc->GetPdf()->plotOn(frame, LineColor(kRed), LineStyle(kDashed), LineWidth(1));
191 var->setVal(-1);
193 list->Add(frame);
194 var->setVal(0);
195 }
196
197 } else {
198 RooCategory *channelCat = (RooCategory *)(&simPdf->indexCat());
199 // TIterator* iter = simPdf->indexCat().typeIterator() ;
200 TIterator *iter = channelCat->typeIterator();
201 RooCatType *tt = NULL;
202 while (nPlots < nPlotsMax && (tt = (RooCatType *)iter->Next())) {
203
204 cout << "on type " << tt->GetName() << " " << endl;
205 // Get pdf associated with state from simpdf
206 RooAbsPdf *pdftmp = simPdf->getPdf(tt->GetName());
207
208 // Generate observables defined by the pdf associated with this state
209 RooArgSet *obstmp = pdftmp->getObservables(*mc->GetObservables());
210 // obstmp->Print();
211
212 obs = ((RooRealVar *)obstmp->first());
213
215 RooRealVar *var = NULL;
216 while (nPlots < nPlotsMax && (var = (RooRealVar *)it.Next())) {
217 TCanvas *c2 = new TCanvas("c2");
218 RooPlot *frame = obs->frame();
219 frame->SetName(Form("frame%d", nPlots));
220 frame->SetYTitle(var->GetName());
221
222 cout << Form("%s==%s::%s", channelCat->GetName(), channelCat->GetName(), tt->GetName()) << endl;
223 cout << tt->GetName() << " " << channelCat->getLabel() << endl;
224 data->plotOn(frame, MarkerSize(1),
225 Cut(Form("%s==%s::%s", channelCat->GetName(), channelCat->GetName(), tt->GetName())),
227
228 Double_t normCount =
229 data->sumEntries(Form("%s==%s::%s", channelCat->GetName(), channelCat->GetName(), tt->GetName()));
230
231 if (strcmp(var->GetName(), "Lumi") == 0) {
232 cout << "working on lumi" << endl;
233 var->setVal(w->var("nominalLumi")->getVal());
234 var->Print();
235 } else {
236 var->setVal(0);
237 }
238 // w->allVars().Print("v");
239 // mc->GetNuisanceParameters()->Print("v");
240 // pdftmp->plotOn(frame,LineWidth(2.));
241 // mc->GetPdf()->plotOn(frame,LineWidth(2.),Slice(*channelCat,tt->GetName()),ProjWData(*data));
242 // pdftmp->plotOn(frame,LineWidth(2.),Slice(*channelCat,tt->GetName()),ProjWData(*data));
243 normCount = pdftmp->expectedEvents(*obs);
244 pdftmp->plotOn(frame, LineWidth(2.), Normalization(normCount, RooAbsReal::NumEvent));
245
246 if (strcmp(var->GetName(), "Lumi") == 0) {
247 cout << "working on lumi" << endl;
248 var->setVal(w->var("nominalLumi")->getVal() + 0.05);
249 var->Print();
250 } else {
251 var->setVal(nSigmaToVary);
252 }
253 // pdftmp->plotOn(frame,LineColor(kRed),LineStyle(kDashed),LineWidth(2));
254 // mc->GetPdf()->plotOn(frame,LineColor(kRed),LineStyle(kDashed),LineWidth(2.),Slice(*channelCat,tt->GetName()),ProjWData(*data));
255 // pdftmp->plotOn(frame,LineColor(kRed),LineStyle(kDashed),LineWidth(2.),Slice(*channelCat,tt->GetName()),ProjWData(*data));
256 normCount = pdftmp->expectedEvents(*obs);
257 pdftmp->plotOn(frame, LineWidth(2.), LineColor(kRed), LineStyle(kDashed),
259
260 if (strcmp(var->GetName(), "Lumi") == 0) {
261 cout << "working on lumi" << endl;
262 var->setVal(w->var("nominalLumi")->getVal() - 0.05);
263 var->Print();
264 } else {
265 var->setVal(-nSigmaToVary);
266 }
267 // pdftmp->plotOn(frame,LineColor(kGreen),LineStyle(kDashed),LineWidth(2));
268 // mc->GetPdf()->plotOn(frame,LineColor(kGreen),LineStyle(kDashed),LineWidth(2),Slice(*channelCat,tt->GetName()),ProjWData(*data));
269 // pdftmp->plotOn(frame,LineColor(kGreen),LineStyle(kDashed),LineWidth(2),Slice(*channelCat,tt->GetName()),ProjWData(*data));
270 normCount = pdftmp->expectedEvents(*obs);
271 pdftmp->plotOn(frame, LineWidth(2.), LineColor(kGreen), LineStyle(kDashed),
273
274 // set them back to normal
275 if (strcmp(var->GetName(), "Lumi") == 0) {
276 cout << "working on lumi" << endl;
277 var->setVal(w->var("nominalLumi")->getVal());
278 var->Print();
279 } else {
280 var->setVal(0);
281 }
282
283 list->Add(frame);
284
285 // quit making plots
286 ++nPlots;
287
288 frame->Draw();
289 c2->SaveAs(Form("%s_%s_%s.pdf", tt->GetName(), obs->GetName(), var->GetName()));
290 delete c2;
291 }
292 }
293 }
294
295 // -------------------------------------------------------
296
297 // now make plots
298 TCanvas *c1 = new TCanvas("c1", "ProfileInspectorDemo", 800, 200);
299 if (list->GetSize() > 4) {
300 double n = list->GetSize();
301 int nx = (int)sqrt(n);
302 int ny = TMath::CeilNint(n / nx);
303 nx = TMath::CeilNint(sqrt(n));
304 c1->Divide(ny, nx);
305 } else
306 c1->Divide(list->GetSize());
307 for (int i = 0; i < list->GetSize(); ++i) {
308 c1->cd(i + 1);
309 list->At(i)->Draw();
310 }
311}
double Double_t
Definition RtypesCore.h:59
@ kRed
Definition Rtypes.h:66
@ kGreen
Definition Rtypes.h:66
@ kDashed
Definition TAttLine.h:48
#define gROOT
Definition TROOT.h:404
char * Form(const char *fmt,...)
R__EXTERN TSystem * gSystem
Definition TSystem.h:559
RooArgSet * getObservables(const RooArgSet &set, Bool_t valueOnly=kTRUE) const
Given a set of possible observables, return the observables that this PDF depends on.
Definition RooAbsArg.h:309
virtual void Print(Option_t *options=0) const
Print the object to the defaultPrintStream().
Definition RooAbsArg.h:337
TIterator * typeIterator() const
const char * getLabel() const
Retrieve current label. Use getCurrentLabel() for more clarity.
RooAbsArg * first() const
virtual void Print(Option_t *options=0) const
This method must be overridden when a class wants to print itself.
TIterator * createIterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:82
virtual const RooArgSet * get() const
Definition RooAbsData.h:128
virtual Double_t sumEntries() const =0
Return effective number of entries in dataset, i.e., sum all weights.
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
virtual RooFitResult * fitTo(RooAbsData &data, 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())
Fit PDF to given dataset.
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:121
virtual Double_t expectedEvents(const RooArgSet *nset) const
Return expected number of events to be used in calculation of extended likelihood.
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:94
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:35
RooCatType is an auxilary class for RooAbsCategory and defines a a single category state.
virtual const Text_t * GetName() const
Returns name of object.
RooCategory is an object to represent discrete states.
Definition RooCategory.h:27
RooDataSet is a container class to hold unbinned data.
Definition RooDataSet.h:36
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:1235
void SetYTitle(const char *title)
Definition RooPlot.cxx:1350
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition RooPlot.cxx:661
RooRealVar represents a variable that can be changed from the outside.
Definition RooRealVar.h:39
virtual void setVal(Double_t value)
Set value of variable to 'value'.
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.
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition ModelConfig.h:30
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return NULL if not existing)
const RooArgSet * GetObservables() const
get RooArgSet for observables (return NULL if not existing)
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
The RooWorkspace is a persistable container for RooFit projects.
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:23
virtual Int_t GetSize() const
Return the capacity of the collection, i.e.
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition TFile.h:54
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:4025
TObject * Next()
Iterator abstract base class.
Definition TIterator.h:30
virtual TObject * Next()=0
A doubly linked list.
Definition TList.h:38
virtual void Add(TObject *obj)
Definition TList.h:81
virtual TObject * At(Int_t idx) const
Returns the object at position idx. Returns 0 if idx is out of range.
Definition TList.cxx:357
virtual const char * GetName() const
Returns name of object.
Definition TNamed.h:47
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition TObject.cxx:200
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition TObject.cxx:267
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:1296
RooCmdArg DataError(Int_t)
RooCmdArg LineWidth(Width_t width)
RooCmdArg MarkerSize(Size_t size)
RooCmdArg Cut(const char *cutSpec)
RooCmdArg Normalization(Double_t scaleFactor)
RooCmdArg LineColor(Color_t color)
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
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition Common.h:18
Namespace for the RooStats classes.
Definition Asimov.h:19
Int_t CeilNint(Double_t x)
Definition TMath.h:649
Definition file.py:1
auto * tt
Definition textangle.C:16