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 // Normally this would be run on the command line
63 cout << "will run standard hist2workspace example" << endl;
64 gROOT->ProcessLine(".! prepareHistFactory .");
65 gROOT->ProcessLine(".! hist2workspace config/example.xml");
66 cout << "\n\n---------------------" << endl;
67 cout << "Done creating example input" << endl;
68 cout << "---------------------\n\n" << endl;
69 }
70
71 } else
72 filename = infile;
73
74 // Try to open the file
75 TFile *file = TFile::Open(filename);
76
77 // if input file was specified byt not found, quit
78 if (!file) {
79 cout << "StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
80 return;
81 }
82
83 // -------------------------------------------------------
84 // Tutorial starts here
85 // -------------------------------------------------------
86
87 // get the workspace out of the file
88 RooWorkspace *w = (RooWorkspace *)file->Get(workspaceName);
89 if (!w) {
90 cout << "workspace not found" << endl;
91 return;
92 }
93
94 // get the modelConfig out of the file
95 ModelConfig *mc = (ModelConfig *)w->obj(modelConfigName);
96
97 // get the modelConfig out of the file
98 RooAbsData *data = w->data(dataName);
99
100 // make sure ingredients are found
101 if (!data || !mc) {
102 w->Print();
103 cout << "data or ModelConfig was not found" << endl;
104 return;
105 }
106
107 // -------------------------------------------------------
108 // create and use the FeldmanCousins tool
109 // to find and plot the 95% confidence interval
110 // on the parameter of interest as specified
111 // in the model config
112 FeldmanCousins fc(*data, *mc);
113 fc.SetConfidenceLevel(0.95); // 95% interval
114 // fc.AdditionalNToysFactor(0.1); // to speed up the result
115 fc.UseAdaptiveSampling(true); // speed it up a bit
116 fc.SetNBins(10); // set how many points per parameter of interest to scan
117 fc.CreateConfBelt(true); // save the information in the belt for plotting
118
119 // Since this tool needs to throw toy MC the PDF needs to be
120 // extended or the tool needs to know how many entries in a dataset
121 // per pseudo experiment.
122 // In the 'number counting form' where the entries in the dataset
123 // are counts, and not values of discriminating variables, the
124 // datasets typically only have one entry and the PDF is not
125 // extended.
126 if (!mc->GetPdf()->canBeExtended()) {
127 if (data->numEntries() == 1)
128 fc.FluctuateNumDataEntries(false);
129 else
130 cout << "Not sure what to do about this model" << endl;
131 }
132
133 // We can use PROOF to speed things along in parallel
134 // ProofConfig pc(*w, 1, "workers=4", kFALSE);
135 // ToyMCSampler* toymcsampler = (ToyMCSampler*) fc.GetTestStatSampler();
136 // toymcsampler->SetProofConfig(&pc); // enable proof
137
138 // Now get the interval
139 PointSetInterval *interval = fc.GetInterval();
140 ConfidenceBelt *belt = fc.GetConfidenceBelt();
141
142 // print out the interval on the first Parameter of Interest
143 RooRealVar *firstPOI = (RooRealVar *)mc->GetParametersOfInterest()->first();
144 cout << "\n95% interval on " << firstPOI->GetName() << " is : [" << interval->LowerLimit(*firstPOI) << ", "
145 << interval->UpperLimit(*firstPOI) << "] " << endl;
146
147 // ---------------------------------------------
148 // No nice plots yet, so plot the belt by hand
149
150 // Ask the calculator which points were scanned
151 RooDataSet *parameterScan = (RooDataSet *)fc.GetPointsToScan();
152 RooArgSet *tmpPoint;
153
154 // make a histogram of parameter vs. threshold
155 TH1F *histOfThresholds =
156 new TH1F("histOfThresholds", "", parameterScan->numEntries(), firstPOI->getMin(), firstPOI->getMax());
157
158 // loop through the points that were tested and ask confidence belt
159 // what the upper/lower thresholds were.
160 // For FeldmanCousins, the lower cut off is always 0
161 for (Int_t i = 0; i < parameterScan->numEntries(); ++i) {
162 tmpPoint = (RooArgSet *)parameterScan->get(i)->clone("temp");
163 double arMax = belt->GetAcceptanceRegionMax(*tmpPoint);
164 double arMin = belt->GetAcceptanceRegionMax(*tmpPoint);
165 double poiVal = tmpPoint->getRealValue(firstPOI->GetName());
166 histOfThresholds->Fill(poiVal, arMax);
167 }
168 histOfThresholds->SetMinimum(0);
169 histOfThresholds->Draw();
170}
int Int_t
Definition RtypesCore.h:45
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
R__EXTERN TSystem * gSystem
Definition TSystem.h:555
double getRealValue(const char *name, double defVal=0.0, bool verbose=false) const
Get value of a RooAbsReal stored in set with given name.
RooAbsArg * first() const
Abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:57
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
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.
virtual double getMin(const char *name=nullptr) const
Get minimum of currently defined range.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
TObject * clone(const char *newname) const override
Definition RooArgSet.h:148
Container class to hold unbinned data.
Definition RooDataSet.h:57
const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
Variable that can be changed from the outside.
Definition RooRealVar.h:37
ConfidenceBelt is a concrete implementation of the ConfInterval interface.
double GetAcceptanceRegionMax(RooArgSet &, double cl=-1., double 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:35
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return nullptr if not existing)
RooAbsPdf * GetPdf() const
get model PDF (return nullptr if pdf has not been specified or does not exist)
PointSetInterval is a concrete implementation of the ConfInterval interface.
double UpperLimit(RooRealVar &param)
return upper limit on a given parameter
double LowerLimit(RooRealVar &param)
return lower limit on a given parameter
Persistable container for RooFit projects.
TObject * Get(const char *namecycle) override
Return pointer to object identified by namecycle.
A ROOT file is composed of a header, followed by consecutive data records (TKey instances) with a wel...
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:4067
1-D histogram with a float per channel (see TH1 documentation)
Definition TH1.h:621
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3340
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition TH1.cxx:3062
virtual void SetMinimum(Double_t minimum=-1111)
Definition TH1.h:404
const char * GetName() const override
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 JSONIO.h:26
Namespace for the RooStats classes.
Definition Asimov.h:19