Logo ROOT   6.21/01
Reference Guide
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 
52 #include "RooStats/ModelConfig.h"
54 
55 using namespace RooFit;
56 using namespace RooStats;
57 using namespace std;
58 
59 void 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);
192  mc->GetPdf()->plotOn(frame, LineColor(kGreen), LineStyle(kDashed), LineWidth(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),
258  Normalization(normCount, RooAbsReal::NumEvent));
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),
272  Normalization(normCount, RooAbsReal::NumEvent));
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 }
virtual Double_t sumEntries() const =0
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
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:1279
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
auto * tt
Definition: textangle.C:16
const RooArgSet * GetObservables() const
get RooArgSet for observables (return NULL if not existing)
Definition: ModelConfig.h:249
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:549
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:3925
RooCmdArg Cut(const char *cutSpec)
void SetName(const char *name)
Set the name of the RooPlot to &#39;name&#39;.
Definition: RooPlot.cxx:1236
RooArgSet * getObservables(const RooArgSet &set, Bool_t valueOnly=kTRUE) const
Return the observables of this pdf given a set of observables.
Definition: RooAbsArg.h:240
RooCmdArg LineColor(Color_t color)
virtual const RooArgSet * get() const
Definition: RooAbsData.h:87
return c1
Definition: legend1.C:41
Definition: Rtypes.h:64
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
void SetYTitle(const char *title)
Definition: RooPlot.cxx:1314
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format...
Definition: TFile.h:48
#define gROOT
Definition: TROOT.h:415
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition: TObject.cxx:195
Definition: Rtypes.h:64
STL namespace.
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
Iterator abstract base class.
Definition: TIterator.h:30
double sqrt(double)
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:128
RooCmdArg MarkerSize(Size_t size)
RooCmdArg DataError(Int_t)
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
RooCmdArg LineStyle(Style_t style)
virtual void Print(Option_t *options=0) const
Print the object to the defaultPrintStream().
Definition: RooAbsArg.h:281
TIterator * createIterator(Bool_t dir=kIterForward) const R__SUGGEST_ALTERNATIVE("begin()
TIterator-style iteration over contained elements.
RooCatType is an auxilary class for RooAbsCategory and defines a a single category state...
Definition: RooCatType.h:22
virtual void Print(Option_t *options=0) const
This method must be overridden when a class wants to print itself.
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:3273
RooRealVar represents a fundamental (non-derived) real-valued object.
Definition: RooRealVar.h:36
A doubly linked list.
Definition: TList.h:44
RooAbsData * data(const char *name) const
Retrieve dataset (binned or unbinned) with given name. A null pointer is returned if not found...
virtual void setVal(Double_t value)
Set value of variable to &#39;value&#39;.
Definition: RooRealVar.cxx:252
R__EXTERN TSystem * gSystem
Definition: TSystem.h:557
const RooAbsCategoryLValue & indexCat() const
RooAbsArg * first() const
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:89
char * Form(const char *fmt,...)
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:44
virtual TObject * At(Int_t idx) const
Returns the object at position idx. Returns 0 if idx is out of range.
Definition: TList.cxx:354
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:31
RooCategory represents a fundamental (non-derived) discrete value object.
Definition: RooCategory.h:24
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, but no plot contents.
TObject * obj(const char *name) const
Return any type of object (RooAbsArg, RooAbsData or generic object) with given name) ...
A RooPlot is a plot frame and a container for graphics objects within that frame. ...
Definition: RooPlot.h:44
The Canvas class.
Definition: TCanvas.h:31
Namespace for the RooStats classes.
Definition: Asimov.h:20
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
Definition: ModelConfig.h:234
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
Definition: ModelConfig.h:237
return c2
Definition: legend2.C:14
double Double_t
Definition: RtypesCore.h:55
RooCmdArg Normalization(Double_t scaleFactor)
RooCmdArg LineWidth(Width_t width)
virtual const char * getLabel() const
Return label string of current state.
Definition: RooCategory.h:39
RooRealVar * var(const char *name) const
Retrieve real-valued variable (RooRealVar) with given name. A null pointer is returned if not found...
RooAbsPdf * getPdf(const char *catName) const
Return the p.d.f associated with the given index category name.
RooAbsPdf, the base class of all PDFs
Definition: RooAbsPdf.h:40
virtual void Add(TObject *obj)
Definition: TList.h:87
Definition: file.py:1
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return NULL if not existing)
Definition: ModelConfig.h:240
virtual TObject * Next()=0
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.
Definition: RooAbsPdf.cxx:1256
virtual Int_t GetSize() const
Return the capacity of the collection, i.e.
Definition: TCollection.h:182
void Print(Option_t *opts=0) const
Print contents of the workspace.
RooSimultaneous facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset...
const Int_t n
Definition: legend1.C:16
Int_t CeilNint(Double_t x)
Definition: TMath.h:689
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:712