Logo ROOT  
Reference Guide
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 
40 #include "RooStats/ModelConfig.h"
42 #include "RooStats/ToyMCSampler.h"
45 
46 using namespace RooFit;
47 using namespace RooStats;
48 
49 void 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 }
RooWorkspace::data
RooAbsData * data(const char *name) const
Retrieve dataset (binned or unbinned) with given name. A null pointer is returned if not found.
Definition: RooWorkspace.cxx:1368
RooWorkspace.h
TH1F.h
PointSetInterval.h
RooAbsCollection::first
RooAbsArg * first() const
Definition: RooAbsCollection.h:177
RooAbsData
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:49
RooAbsRealLValue::getMax
virtual Double_t getMax(const char *name=0) const
Get maximum of currently defined range.
Definition: RooAbsRealLValue.h:89
RooStats::PointSetInterval::UpperLimit
Double_t UpperLimit(RooRealVar &param)
return upper limit on a given parameter
Definition: PointSetInterval.cxx:147
RooArgSet::getRealValue
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.
Definition: RooArgSet.cxx:343
TH1::SetMinimum
virtual void SetMinimum(Double_t minimum=-1111)
Definition: TH1.h:399
fc
static struct mg_connection * fc(struct mg_context *ctx)
Definition: civetweb.c:3728
RooArgSet::clone
TObject * clone(const char *newname) const override
Definition: RooArgSet.h:83
ConfidenceBelt.h
RooStats::FeldmanCousins
The FeldmanCousins class (like the Feldman-Cousins technique) is essentially a specific configuration...
Definition: FeldmanCousins.h:33
TFile::Open
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
ToyMCSampler.h
TSystem::AccessPathName
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
TFile.h
RooAbsRealLValue::getMin
virtual Double_t getMin(const char *name=0) const
Get miniminum of currently defined range.
Definition: RooAbsRealLValue.h:86
RooDataSet::get
virtual const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
Definition: RooDataSet.cxx:1051
TROOT.h
TSystem.h
RooWorkspace::Print
void Print(Option_t *opts=0) const
Print contents of the workspace.
Definition: RooWorkspace.cxx:2194
ModelConfig.h
RooFit
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition: RooCFunction1Binding.h:29
RooStats::ConfidenceBelt
ConfidenceBelt is a concrete implementation of the ConfInterval interface.
Definition: ConfidenceBelt.h:156
TH1::Fill
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3346
RooWorkspace::obj
TObject * obj(const char *name) const
Return any type of object (RooAbsArg, RooAbsData or generic object) with given name)
Definition: RooWorkspace.cxx:2106
RooAbsData::numEntries
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Definition: RooAbsData.cxx:307
TFile
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:54
RooStats::ModelConfig::GetPdf
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
Definition: ModelConfig.h:234
gSystem
R__EXTERN TSystem * gSystem
Definition: TSystem.h:559
FeldmanCousins.h
RooStats::ModelConfig::GetParametersOfInterest
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
Definition: ModelConfig.h:237
RooWorkspace
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
RooAbsData.h
RooStats
Namespace for the RooStats classes.
Definition: Asimov.h:19
file
Definition: file.py:1
TH1F
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:575
RooStats::ConfidenceBelt::GetAcceptanceRegionMax
Double_t GetAcceptanceRegionMax(RooArgSet &, Double_t cl=-1., Double_t leftside=-1.)
Definition: ConfidenceBelt.cxx:101
RooDataSet
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
RooRealVar
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:39
RooStats::PointSetInterval
PointSetInterval is a concrete implementation of the ConfInterval interface.
Definition: PointSetInterval.h:21
RooStats::ModelConfig
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
RooStats::PointSetInterval::LowerLimit
Double_t LowerLimit(RooRealVar &param)
return lower limit on a given parameter
Definition: PointSetInterval.cxx:160
RooAbsPdf::canBeExtended
Bool_t canBeExtended() const
If true, PDF can provide extended likelihood term.
Definition: RooAbsPdf.h:238
RooArgSet
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:29
gROOT
#define gROOT
Definition: TROOT.h:406
int
TH1::Draw
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:3069