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 // Normally this would be run on the command line
80 cout << "will run standard hist2workspace example" << endl;
81 gROOT->ProcessLine(".! prepareHistFactory .");
82 gROOT->ProcessLine(".! hist2workspace config/example.xml");
83 cout << "\n\n---------------------" << endl;
84 cout << "Done creating example input" << endl;
85 cout << "---------------------\n\n" << endl;
86 }
87
88 } else
89 filename = infile;
90
91 // Try to open the file
92 TFile *file = TFile::Open(filename);
93
94 // -------------------------------------------------------
95 // Tutorial starts here
96 // -------------------------------------------------------
97
98 // get the workspace out of the file
99 RooWorkspace *w = (RooWorkspace *)file->Get(workspaceName);
100
101 // get the modelConfig out of the file
102 ModelConfig *mc = (ModelConfig *)w->obj(modelConfigName);
103
104 // get the modelConfig out of the file
105 RooAbsData *data = w->data(dataName);
106
107 // ---------------------------------------------
108 // create and use the ProfileLikelihoodCalculator
109 // to find and plot the 95% confidence interval
110 // on the parameter of interest as specified
111 // in the model config
113 pl.SetConfidenceLevel(confLevel); // 95% interval
114 LikelihoodInterval *interval = pl.GetInterval();
115
116 // print out the interval on the first Parameter of Interest
117 RooRealVar *firstPOI = (RooRealVar *)mc->GetParametersOfInterest()->first();
118 cout << "\n>>>> RESULT : " << confLevel * 100 << "% interval on " << firstPOI->GetName() << " is : ["
119 << interval->LowerLimit(*firstPOI) << ", " << interval->UpperLimit(*firstPOI) << "]\n " << endl;
120
121 // make a plot
122
123 cout << "making a plot of the profile likelihood function ....(if it is taking a lot of time use less points or the "
124 "TF1 drawing option)\n";
126 plot.SetNPoints(nScanPoints); // do not use too many points, it could become very slow for some models
127 if (poiXMin < poiXMax)
128 plot.SetRange(poiXMin, poiXMax);
129 TString opt;
130 if (plotAsTF1)
131 opt += TString("tf1");
132 plot.Draw(opt); // use option TF1 if too slow (plot.Draw("tf1")
133
134 // if requested perform also an hypothesis test for the significance
135 if (doHypoTest) {
136 RooArgSet nullparams("nullparams");
137 nullparams.addClone(*firstPOI);
138 nullparams.setRealValue(firstPOI->GetName(), nullParamValue);
139 pl.SetNullParameters(nullparams);
140 std::cout << "Perform Test of Hypothesis : null Hypothesis is " << firstPOI->GetName() << nullParamValue
141 << std::endl;
142 auto result = pl.GetHypoTest();
143 std::cout << "\n>>>> Hypotheis Test Result \n";
144 result->Print();
145 }
146}
winID h TVirtualViewer3D TVirtualGLPainter char TVirtualGLPainter plot
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char filename
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
#define gROOT
Definition TROOT.h:407
R__EXTERN TSystem * gSystem
Definition TSystem.h:555
RooAbsArg * first() const
Abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:57
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
Variable that can be changed from the outside.
Definition RooRealVar.h:37
This class provides simple and straightforward utilities to plot a LikelihoodInterval object.
LikelihoodInterval is a concrete implementation of the RooStats::ConfInterval interface.
double UpperLimit(const RooRealVar &param)
return the upper bound of the interval on a given parameter
double LowerLimit(const RooRealVar &param)
return the lower 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:35
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return nullptr if not existing)
The ProfileLikelihoodCalculator is a concrete implementation of CombinedCalculator (the interface cla...
Persistable container for RooFit projects.
TObject * Get(const char *namecycle) override
Return pointer to object identified by namecycle.
A ROOT file is an on-disk file, usually with extension .root, that stores objects in a file-system-li...
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:4082
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition TObject.cxx:274
Basic string class.
Definition TString.h:139
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
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition JSONIO.h:26
Namespace for the RooStats classes.
Definition Asimov.h:19