Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
41
42using namespace RooFit;
43using namespace RooStats;
44
45struct 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
56ProfileLikelihoodOptions optPL;
57
58void 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}
#define gROOT
Definition TROOT.h:406
R__EXTERN TSystem * gSystem
Definition TSystem.h:559
RooAbsArg * first() const
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:49
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:29
RooRealVar represents a variable that can be changed from the outside.
Definition RooRealVar.h:39
This class provides simple and straightforward utilities to plot a LikelihoodInterval object.
LikelihoodInterval is a concrete implementation of the RooStats::ConfInterval interface.
Double_t LowerLimit(const RooRealVar &param)
return the lower bound of the interval on a given parameter
Double_t UpperLimit(const RooRealVar &param)
return the upper bound of the interval on a given parameter
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)
The ProfileLikelihoodCalculator is a concrete implementation of CombinedCalculator (the interface cla...
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.
TObject * obj(const char *name) const
Return any type of object (RooAbsArg, RooAbsData or generic object) with given name)
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:3997
virtual const char * GetName() const
Returns name of object.
Definition TNamed.h:47
Basic string class.
Definition TString.h:136
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:1294
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
Definition file.py:1