Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
StandardTestStatDistributionDemo.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_roostats
3/// \notebook
4/// StandardTestStatDistributionDemo.C
5///
6/// This simple script plots the sampling distribution of the profile likelihood
7/// ratio test statistic based on the input Model File. To do this one needs to
8/// specify the value of the parameter of interest that will be used for evaluating
9/// the test statistic and the value of the parameters used for generating the toy data.
10/// In this case, it uses the upper-limit estimated from the ProfileLikleihoodCalculator,
11/// which assumes the asymptotic chi-square distribution for -2 log profile likelihood ratio.
12/// Thus, the script is handy for checking to see if the asymptotic approximations are valid.
13/// To aid, that comparison, the script overlays a chi-square distribution as well.
14/// The most common parameter of interest is a parameter proportional to the signal rate,
15/// and often that has a lower-limit of 0, which breaks the standard chi-square distribution.
16/// Thus the script allows the parameter to be negative so that the overlay chi-square is
17/// the correct asymptotic distribution.
18///
19/// \macro_image
20/// \macro_output
21/// \macro_code
22///
23/// \author Kyle Cranmer
24
25#include "TFile.h"
26#include "TROOT.h"
27#include "TH1F.h"
28#include "TCanvas.h"
29#include "TSystem.h"
30#include "TF1.h"
31#include "TSystem.h"
32
33#include "RooWorkspace.h"
34#include "RooAbsData.h"
35
41
47
48using namespace RooFit;
49using namespace RooStats;
50
51bool useProof = false; // flag to control whether to use Proof
52int nworkers = 0; // number of workers (default use all available cores)
53
54// -------------------------------------------------------
55// The actual macro
56
57void StandardTestStatDistributionDemo(const char *infile = "", const char *workspaceName = "combined",
58 const char *modelConfigName = "ModelConfig", const char *dataName = "obsData")
59{
60
61 // the number of toy MC used to generate the distribution
62 int nToyMC = 1000;
63 // The parameter below is needed for asymptotic distribution to be chi-square,
64 // but set to false if your model is not numerically stable if mu<0
65 bool allowNegativeMu = true;
66
67 // -------------------------------------------------------
68 // First part is just to access a user-defined file
69 // or create the standard example file if it doesn't exist
70 const char *filename = "";
71 if (!strcmp(infile, "")) {
72 filename = "results/example_combined_GaussExample_model.root";
73 bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code
74 // if file does not exists generate with histfactory
75 if (!fileExist) {
76 // Normally this would be run on the command line
77 cout << "will run standard hist2workspace example" << endl;
78 gROOT->ProcessLine(".! prepareHistFactory .");
79 gROOT->ProcessLine(".! hist2workspace config/example.xml");
80 cout << "\n\n---------------------" << endl;
81 cout << "Done creating example input" << endl;
82 cout << "---------------------\n\n" << endl;
83 }
84
85 } else
86 filename = infile;
87
88 // Try to open the file
89 TFile *file = TFile::Open(filename);
90
91 // if input file was specified byt not found, quit
92 if (!file) {
93 cout << "StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
94 return;
95 }
96
97 // -------------------------------------------------------
98 // Now get the data and workspace
99
100 // get the workspace out of the file
101 RooWorkspace *w = (RooWorkspace *)file->Get(workspaceName);
102 if (!w) {
103 cout << "workspace not found" << endl;
104 return;
105 }
106
107 // get the modelConfig out of the file
108 ModelConfig *mc = (ModelConfig *)w->obj(modelConfigName);
109
110 // get the modelConfig out of the file
111 RooAbsData *data = w->data(dataName);
112
113 // make sure ingredients are found
114 if (!data || !mc) {
115 w->Print();
116 cout << "data or ModelConfig was not found" << endl;
117 return;
118 }
119
120 mc->Print();
121 // -------------------------------------------------------
122 // Now find the upper limit based on the asymptotic results
123 RooRealVar *firstPOI = (RooRealVar *)mc->GetParametersOfInterest()->first();
125 LikelihoodInterval *interval = plc.GetInterval();
126 double plcUpperLimit = interval->UpperLimit(*firstPOI);
127 delete interval;
128 cout << "\n\n--------------------------------------" << endl;
129 cout << "Will generate sampling distribution at " << firstPOI->GetName() << " = " << plcUpperLimit << endl;
130 int nPOI = mc->GetParametersOfInterest()->getSize();
131 if (nPOI > 1) {
132 cout << "not sure what to do with other parameters of interest, but here are their values" << endl;
133 mc->GetParametersOfInterest()->Print("v");
134 }
135
136 // -------------------------------------------------------
137 // create the test stat sampler
139
140 // to avoid effects from boundary and simplify asymptotic comparison, set min=-max
141 if (allowNegativeMu)
142 firstPOI->setMin(-1 * firstPOI->getMax());
143
144 // temporary RooArgSet
145 RooArgSet poi;
146 poi.add(*mc->GetParametersOfInterest());
147
148 // create and configure the ToyMCSampler
149 ToyMCSampler sampler(ts, nToyMC);
150 sampler.SetPdf(*mc->GetPdf());
151 sampler.SetObservables(*mc->GetObservables());
152 sampler.SetGlobalObservables(*mc->GetGlobalObservables());
153 if (!mc->GetPdf()->canBeExtended() && (data->numEntries() == 1)) {
154 cout << "tell it to use 1 event" << endl;
155 sampler.SetNEventsPerToy(1);
156 }
157 firstPOI->setVal(plcUpperLimit); // set POI value for generation
158 sampler.SetParametersForTestStat(*mc->GetParametersOfInterest()); // set POI value for evaluation
159
160 if (useProof) {
161 ProofConfig pc(*w, nworkers, "", false);
162 sampler.SetProofConfig(&pc); // enable proof
163 }
164
165 firstPOI->setVal(plcUpperLimit);
166 RooArgSet allParameters;
167 allParameters.add(*mc->GetParametersOfInterest());
168 allParameters.add(*mc->GetNuisanceParameters());
169 allParameters.Print("v");
170
171 SamplingDistribution *sampDist = sampler.GetSamplingDistribution(allParameters);
173 plot.AddSamplingDistribution(sampDist);
174 plot.GetTH1F(sampDist)->GetYaxis()->SetTitle(
175 Form("f(-log #lambda(#mu=%.2f) | #mu=%.2f)", plcUpperLimit, plcUpperLimit));
176 plot.SetAxisTitle(Form("-log #lambda(#mu=%.2f)", plcUpperLimit));
177
178 TCanvas *c1 = new TCanvas("c1");
179 c1->SetLogy();
180 plot.Draw();
181 double min = plot.GetTH1F(sampDist)->GetXaxis()->GetXmin();
182 double max = plot.GetTH1F(sampDist)->GetXaxis()->GetXmax();
183
184 TF1 *f = new TF1("f", Form("2*ROOT::Math::chisquared_pdf(2*x,%d,0)", nPOI), min, max);
185 f->Draw("same");
186 c1->SaveAs("standard_test_stat_distribution.pdf");
187}
#define f(i)
Definition RSha256.hxx:104
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
#define gROOT
Definition TROOT.h:406
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition TString.cxx:2489
R__EXTERN TSystem * gSystem
Definition TSystem.h:555
Int_t getSize() const
Return the number of elements in the collection.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
RooAbsArg * first() const
void Print(Option_t *options=nullptr) const override
This method must be overridden when a class wants to print itself.
Abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:57
bool canBeExtended() const
If true, PDF can provide extended likelihood term.
Definition RooAbsPdf.h:219
virtual double getMax(const char *name=nullptr) const
Get maximum of currently defined range.
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
void setVal(double value) override
Set value of variable to 'value'.
void setMin(const char *name, double value)
Set minimum of name range to given value.
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
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition ModelConfig.h:35
const RooArgSet * GetGlobalObservables() const
get RooArgSet for global observables (return nullptr if not existing)
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return nullptr if not existing)
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return nullptr if not existing)
void Print(Option_t *option="") const override
overload the print method
const RooArgSet * GetObservables() const
get RooArgSet for observables (return nullptr if not existing)
RooAbsPdf * GetPdf() const
get model PDF (return nullptr if pdf has not been specified or does not exist)
The ProfileLikelihoodCalculator is a concrete implementation of CombinedCalculator (the interface cla...
ProfileLikelihoodTestStat is an implementation of the TestStatistic interface that calculates the pro...
Holds configuration options for proof and proof-lite.
Definition ProofConfig.h:45
This class provides simple and straightforward utilities to plot SamplingDistribution objects.
This class simply holds a sampling distribution of some test statistic.
const std::vector< double > & GetSamplingDistribution() const
Get test statistics values.
ToyMCSampler is an implementation of the TestStatSampler interface.
Persistable container for RooFit projects.
The Canvas class.
Definition TCanvas.h:23
TObject * Get(const char *namecycle) override
Return pointer to object identified by namecycle.
1-Dim function class
Definition TF1.h:233
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:4089
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
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
return c1
Definition legend1.C:41
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
__device__ AFloat max(AFloat x, AFloat y)
Definition Kernels.cuh:207