Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
StandardFeldmanCousinsDemo.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_roostats
3/// \notebook -js
4/// Standard demo of the Feldman-Cousins calculator
5/// StandardFeldmanCousinsDemo
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 FeldmanCousins tools is a classical frequentist calculation
21/// based on the Neyman Construction. The test statistic can be
22/// generalized for nuisance parameters by using the profile likelihood ratio.
23/// But unlike the ProfileLikelihoodCalculator, this tool explicitly
24/// builds the sampling distribution of the test statistic via toy Monte Carlo.
25///
26/// \macro_image
27/// \macro_output
28/// \macro_code
29///
30/// \author Kyle Cranmer
31
32#include "TFile.h"
33#include "TROOT.h"
34#include "TH1F.h"
35#include "TSystem.h"
36
37#include "RooWorkspace.h"
38#include "RooAbsData.h"
39
45
46using namespace RooFit;
47using namespace RooStats;
48
49void StandardFeldmanCousinsDemo(const char *infile = "", const char *workspaceName = "combined",
50 const char *modelConfigName = "ModelConfig", const char *dataName = "obsData")
51{
52
53 // -------------------------------------------------------
54 // First part is just to access a user-defined file
55 // or create the standard example file if it doesn't exist
56 const char *filename = "";
57 if (!strcmp(infile, "")) {
58 filename = "results/example_combined_GaussExample_model.root";
59 bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code
60 // if file does not exists generate with histfactory
61 if (!fileExist) {
62#ifdef _WIN32
63 cout << "HistFactory file cannot be generated on Windows - exit" << endl;
64 return;
65#endif
66 // Normally this would be run on the command line
67 cout << "will run standard hist2workspace example" << endl;
68 gROOT->ProcessLine(".! prepareHistFactory .");
69 gROOT->ProcessLine(".! hist2workspace config/example.xml");
70 cout << "\n\n---------------------" << endl;
71 cout << "Done creating example input" << endl;
72 cout << "---------------------\n\n" << endl;
73 }
74
75 } else
76 filename = infile;
77
78 // Try to open the file
79 TFile *file = TFile::Open(filename);
80
81 // if input file was specified byt not found, quit
82 if (!file) {
83 cout << "StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
84 return;
85 }
86
87 // -------------------------------------------------------
88 // Tutorial starts here
89 // -------------------------------------------------------
90
91 // get the workspace out of the file
92 RooWorkspace *w = (RooWorkspace *)file->Get(workspaceName);
93 if (!w) {
94 cout << "workspace not found" << endl;
95 return;
96 }
97
98 // get the modelConfig out of the file
99 ModelConfig *mc = (ModelConfig *)w->obj(modelConfigName);
100
101 // get the modelConfig out of the file
102 RooAbsData *data = w->data(dataName);
103
104 // make sure ingredients are found
105 if (!data || !mc) {
106 w->Print();
107 cout << "data or ModelConfig was not found" << endl;
108 return;
109 }
110
111 // -------------------------------------------------------
112 // create and use the FeldmanCousins tool
113 // to find and plot the 95% confidence interval
114 // on the parameter of interest as specified
115 // in the model config
116 FeldmanCousins fc(*data, *mc);
117 fc.SetConfidenceLevel(0.95); // 95% interval
118 // fc.AdditionalNToysFactor(0.1); // to speed up the result
119 fc.UseAdaptiveSampling(true); // speed it up a bit
120 fc.SetNBins(10); // set how many points per parameter of interest to scan
121 fc.CreateConfBelt(true); // save the information in the belt for plotting
122
123 // Since this tool needs to throw toy MC the PDF needs to be
124 // extended or the tool needs to know how many entries in a dataset
125 // per pseudo experiment.
126 // In the 'number counting form' where the entries in the dataset
127 // are counts, and not values of discriminating variables, the
128 // datasets typically only have one entry and the PDF is not
129 // extended.
130 if (!mc->GetPdf()->canBeExtended()) {
131 if (data->numEntries() == 1)
132 fc.FluctuateNumDataEntries(false);
133 else
134 cout << "Not sure what to do about this model" << endl;
135 }
136
137 // We can use PROOF to speed things along in parallel
138 // ProofConfig pc(*w, 1, "workers=4", kFALSE);
139 // ToyMCSampler* toymcsampler = (ToyMCSampler*) fc.GetTestStatSampler();
140 // toymcsampler->SetProofConfig(&pc); // enable proof
141
142 // Now get the interval
143 PointSetInterval *interval = fc.GetInterval();
144 ConfidenceBelt *belt = fc.GetConfidenceBelt();
145
146 // print out the interval on the first Parameter of Interest
147 RooRealVar *firstPOI = (RooRealVar *)mc->GetParametersOfInterest()->first();
148 cout << "\n95% interval on " << firstPOI->GetName() << " is : [" << interval->LowerLimit(*firstPOI) << ", "
149 << interval->UpperLimit(*firstPOI) << "] " << endl;
150
151 // ---------------------------------------------
152 // No nice plots yet, so plot the belt by hand
153
154 // Ask the calculator which points were scanned
155 RooDataSet *parameterScan = (RooDataSet *)fc.GetPointsToScan();
156 RooArgSet *tmpPoint;
157
158 // make a histogram of parameter vs. threshold
159 TH1F *histOfThresholds =
160 new TH1F("histOfThresholds", "", parameterScan->numEntries(), firstPOI->getMin(), firstPOI->getMax());
161
162 // loop through the points that were tested and ask confidence belt
163 // what the upper/lower thresholds were.
164 // For FeldmanCousins, the lower cut off is always 0
165 for (Int_t i = 0; i < parameterScan->numEntries(); ++i) {
166 tmpPoint = (RooArgSet *)parameterScan->get(i)->clone("temp");
167 double arMax = belt->GetAcceptanceRegionMax(*tmpPoint);
168 double arMin = belt->GetAcceptanceRegionMax(*tmpPoint);
169 double poiVal = tmpPoint->getRealValue(firstPOI->GetName());
170 histOfThresholds->Fill(poiVal, arMax);
171 }
172 histOfThresholds->SetMinimum(0);
173 histOfThresholds->Draw();
174}
int Int_t
Definition RtypesCore.h:45
#define gROOT
Definition TROOT.h:404
R__EXTERN TSystem * gSystem
Definition TSystem.h:559
static struct mg_connection * fc(struct mg_context *ctx)
Definition civetweb.c:3728
Double_t getRealValue(const char *name, Double_t defVal=0, Bool_t verbose=kFALSE) const
Get value of a RooAbsReal stored in set with given name.
RooAbsArg * first() const
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:82
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Bool_t canBeExtended() const
If true, PDF can provide extended likelihood term.
Definition RooAbsPdf.h:262
virtual Double_t getMax(const char *name=0) const
Get maximum of currently defined range.
virtual Double_t getMin(const char *name=0) const
Get miniminum of currently defined range.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:35
TObject * clone(const char *newname) const override
Definition RooArgSet.h:128
RooDataSet is a container class to hold unbinned data.
Definition RooDataSet.h:36
virtual const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
RooRealVar represents a variable that can be changed from the outside.
Definition RooRealVar.h:39
ConfidenceBelt is a concrete implementation of the ConfInterval interface.
Double_t GetAcceptanceRegionMax(RooArgSet &, Double_t cl=-1., Double_t leftside=-1.)
The FeldmanCousins class (like the Feldman-Cousins technique) is essentially a specific configuration...
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)
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
PointSetInterval is a concrete implementation of the ConfInterval interface.
Double_t UpperLimit(RooRealVar &param)
return upper limit on a given parameter
Double_t LowerLimit(RooRealVar &param)
return lower limit on a given parameter
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:4025
1-D histogram with a float per channel (see TH1 documentation)}
Definition TH1.h:575
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3351
virtual void SetMinimum(Double_t minimum=-1111)
Definition TH1.h:399
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition TH1.cxx:3074
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:1296
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition Common.h:18
Namespace for the RooStats classes.
Definition Asimov.h:19
Definition file.py:1