Logo ROOT   6.21/01
Reference Guide
StandardProfileLikelihoodDemo.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_roostats
3 /// \notebook
4 /// Standard demo of the Profile Likelihood calculator
5 /// StandardProfileLikelihoodDemo
6 ///
7 /// This is a standard demo that can be used with any ROOT file
8 /// prepared in the standard way. You specify:
9 /// - name for input ROOT file
10 /// - name of workspace inside ROOT file that holds model and data
11 /// - name of ModelConfig that specifies details for calculator tools
12 /// - name of dataset
13 ///
14 /// With default parameters the macro will attempt to run the
15 /// standard hist2workspace example and read the ROOT file
16 /// that it produces.
17 ///
18 /// The actual heart of the demo is only about 10 lines long.
19 ///
20 /// The ProfileLikelihoodCalculator is based on Wilks's theorem
21 /// and the asymptotic properties of the profile likelihood ratio
22 /// (eg. that it is chi-square distributed for the true value).
23 ///
24 /// \macro_image
25 /// \macro_output
26 /// \macro_code
27 ///
28 /// \author Kyle Cranmer
29 
30 #include "TFile.h"
31 #include "TROOT.h"
32 #include "TSystem.h"
33 #include "RooWorkspace.h"
34 #include "RooAbsData.h"
35 #include "RooRealVar.h"
36 
37 #include "RooStats/ModelConfig.h"
41 
42 using namespace RooFit;
43 using namespace RooStats;
44 
45 struct ProfileLikelihoodOptions {
46 
47  double confLevel = 0.95;
48  int nScanPoints = 50;
49  bool plotAsTF1 = false;
50  double poiMinPlot = 1;
51  double poiMaxPlot = 0;
52  bool doHypoTest = false;
53  double nullValue = 0;
54 };
55 
56 ProfileLikelihoodOptions optPL;
57 
58 void StandardProfileLikelihoodDemo(const char *infile = "", const char *workspaceName = "combined",
59  const char *modelConfigName = "ModelConfig", const char *dataName = "obsData")
60 {
61 
62  double confLevel = optPL.confLevel;
63  double nScanPoints = optPL.nScanPoints;
64  bool plotAsTF1 = optPL.plotAsTF1;
65  double poiXMin = optPL.poiMinPlot;
66  double poiXMax = optPL.poiMaxPlot;
67  bool doHypoTest = optPL.doHypoTest;
68  double nullParamValue = optPL.nullValue;
69 
70  // -------------------------------------------------------
71  // First part is just to access a user-defined file
72  // or create the standard example file if it doesn't exist
73  const char *filename = "";
74  if (!strcmp(infile, "")) {
75  filename = "results/example_combined_GaussExample_model.root";
76  bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code
77  // if file does not exists generate with histfactory
78  if (!fileExist) {
79 #ifdef _WIN32
80  cout << "HistFactory file cannot be generated on Windows - exit" << endl;
81  return;
82 #endif
83  // Normally this would be run on the command line
84  cout << "will run standard hist2workspace example" << endl;
85  gROOT->ProcessLine(".! prepareHistFactory .");
86  gROOT->ProcessLine(".! hist2workspace config/example.xml");
87  cout << "\n\n---------------------" << endl;
88  cout << "Done creating example input" << endl;
89  cout << "---------------------\n\n" << endl;
90  }
91 
92  } else
93  filename = infile;
94 
95  // Try to open the file
96  TFile *file = TFile::Open(filename);
97 
98  // if input file was specified byt not found, quit
99  if (!file) {
100  cout << "StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
101  return;
102  }
103 
104  // -------------------------------------------------------
105  // Tutorial starts here
106  // -------------------------------------------------------
107 
108  // get the workspace out of the file
109  RooWorkspace *w = (RooWorkspace *)file->Get(workspaceName);
110  if (!w) {
111  cout << "workspace not found" << endl;
112  return;
113  }
114 
115  // get the modelConfig out of the file
116  ModelConfig *mc = (ModelConfig *)w->obj(modelConfigName);
117 
118  // get the modelConfig out of the file
119  RooAbsData *data = w->data(dataName);
120 
121  // make sure ingredients are found
122  if (!data || !mc) {
123  w->Print();
124  cout << "data or ModelConfig was not found" << endl;
125  return;
126  }
127 
128  // ---------------------------------------------
129  // create and use the ProfileLikelihoodCalculator
130  // to find and plot the 95% confidence interval
131  // on the parameter of interest as specified
132  // in the model config
133  ProfileLikelihoodCalculator pl(*data, *mc);
134  pl.SetConfidenceLevel(confLevel); // 95% interval
135  LikelihoodInterval *interval = pl.GetInterval();
136 
137  // print out the interval on the first Parameter of Interest
138  RooRealVar *firstPOI = (RooRealVar *)mc->GetParametersOfInterest()->first();
139  cout << "\n>>>> RESULT : " << confLevel * 100 << "% interval on " << firstPOI->GetName() << " is : ["
140  << interval->LowerLimit(*firstPOI) << ", " << interval->UpperLimit(*firstPOI) << "]\n " << endl;
141 
142  // make a plot
143 
144  cout << "making a plot of the profile likelihood function ....(if it is taking a lot of time use less points or the "
145  "TF1 drawing option)\n";
146  LikelihoodIntervalPlot plot(interval);
147  plot.SetNPoints(nScanPoints); // do not use too many points, it could become very slow for some models
148  if (poiXMin < poiXMax)
149  plot.SetRange(poiXMin, poiXMax);
150  TString opt;
151  if (plotAsTF1)
152  opt += TString("tf1");
153  plot.Draw(opt); // use option TF1 if too slow (plot.Draw("tf1")
154 
155  // if requested perform also an hypothesis test for the significance
156  if (doHypoTest) {
157  RooArgSet nullparams("nullparams");
158  nullparams.addClone(*firstPOI);
159  nullparams.setRealValue(firstPOI->GetName(), nullParamValue);
160  pl.SetNullParameters(nullparams);
161  std::cout << "Perform Test of Hypothesis : null Hypothesis is " << firstPOI->GetName() << nullParamValue
162  << std::endl;
163  auto result = pl.GetHypoTest();
164  std::cout << "\n>>>> Hypotheis Test Result \n";
165  result->Print();
166  }
167 }
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
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
LikelihoodInterval is a concrete implementation of the RooStats::ConfInterval interface.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format...
Definition: TFile.h:48
This class provides simple and straightforward utilities to plot a LikelihoodInterval object...
The ProfileLikelihoodCalculator is a concrete implementation of CombinedCalculator (the interface cla...
#define gROOT
Definition: TROOT.h:415
Basic string class.
Definition: TString.h:131
Double_t LowerLimit(const RooRealVar &param)
return the lower bound of the interval on a given parameter
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
RooRealVar represents a fundamental (non-derived) real-valued object.
Definition: RooRealVar.h:36
RooAbsData * data(const char *name) const
Retrieve dataset (binned or unbinned) with given name. A null pointer is returned if not found...
Double_t UpperLimit(const RooRealVar &param)
return the upper bound of the interval on a given parameter
R__EXTERN TSystem * gSystem
Definition: TSystem.h:557
RooAbsArg * first() const
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:44
TObject * obj(const char *name) const
Return any type of object (RooAbsArg, RooAbsData or generic object) with given name) ...
Namespace for the RooStats classes.
Definition: Asimov.h:20
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
Definition: ModelConfig.h:237
Definition: file.py:1
void Print(Option_t *opts=0) const
Print contents of the workspace.
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43