Logo ROOT  
Reference Guide
StandardFrequentistDiscovery.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_roostats
3/// \notebook
4/// StandardFrequentistDiscovery
5///
6/// This is a standard demo that can be used with any ROOT file
7/// prepared in the standard way. You specify:
8/// - name for input ROOT file
9/// - name of workspace inside ROOT file that holds model and data
10/// - name of ModelConfig that specifies details for calculator tools
11/// - name of dataset
12///
13/// With default parameters the macro will attempt to run the
14/// standard hist2workspace example and read the ROOT file
15/// that it produces.
16///
17/// \macro_image
18/// \macro_output
19/// \macro_code
20///
21/// \authors Sven Kreiss, Kyle Cranmer
22
23#include "TFile.h"
24#include "TROOT.h"
25#include "TH1F.h"
26#include "TF1.h"
27#include "TCanvas.h"
28#include "TStopwatch.h"
29
30#include "RooWorkspace.h"
31#include "RooAbsData.h"
32#include "RooRandom.h"
33#include "RooRealSumPdf.h"
34#include "RooNumIntConfig.h"
35
46
48#include "TSystem.h"
49
50#include <vector>
51
52using namespace RooFit;
53using namespace RooStats;
54
55double StandardFrequentistDiscovery(const char *infile = "", const char *workspaceName = "channel1",
56 const char *modelConfigNameSB = "ModelConfig", const char *dataName = "obsData",
57 int toys = 1000, double poiValueForBackground = 0.0, double poiValueForSignal = 1.0)
58{
59
60 // The workspace contains the model for s+b. The b model is "autogenerated"
61 // by copying s+b and setting the one parameter of interest to zero.
62 // To keep the script simple, multiple parameters of interest or different
63 // functional forms of the b model are not supported.
64
65 // for now, assume there is only one parameter of interest, and these are
66 // its values:
67
68 // -------------------------------------------------------
69 // First part is just to access a user-defined file
70 // or create the standard example file if it doesn't exist
71 const char *filename = "";
72 if (!strcmp(infile, "")) {
73 filename = "results/example_channel1_GammaExample_model.root";
74 bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code
75 // if file does not exists generate with histfactory
76 if (!fileExist) {
77#ifdef _WIN32
78 cout << "HistFactory file cannot be generated on Windows - exit" << endl;
79 return -1;
80#endif
81 // Normally this would be run on the command line
82 cout << "will run standard hist2workspace example" << endl;
83 gROOT->ProcessLine(".! prepareHistFactory .");
84 gROOT->ProcessLine(".! hist2workspace config/example.xml");
85 cout << "\n\n---------------------" << endl;
86 cout << "Done creating example input" << endl;
87 cout << "---------------------\n\n" << endl;
88 }
89
90 } else
91 filename = infile;
92
93 // Try to open the file
94 TFile *file = TFile::Open(filename);
95
96 // if input file was specified byt not found, quit
97 if (!file) {
98 cout << "StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
99 return -1;
100 }
101
102 // -------------------------------------------------------
103 // Tutorial starts here
104 // -------------------------------------------------------
105
106 TStopwatch *mn_t = new TStopwatch;
107 mn_t->Start();
108
109 // get the workspace out of the file
110 RooWorkspace *w = (RooWorkspace *)file->Get(workspaceName);
111 if (!w) {
112 cout << "workspace not found" << endl;
113 return -1.0;
114 }
115
116 // get the modelConfig out of the file
117 ModelConfig *mc = (ModelConfig *)w->obj(modelConfigNameSB);
118
119 // get the data out of the file
120 RooAbsData *data = w->data(dataName);
121
122 // make sure ingredients are found
123 if (!data || !mc) {
124 w->Print();
125 cout << "data or ModelConfig was not found" << endl;
126 return -1.0;
127 }
128
129 RooRealVar *firstPOI = (RooRealVar *)mc->GetParametersOfInterest()->first();
130 firstPOI->setVal(poiValueForSignal);
132 // create null model
133 ModelConfig *mcNull = mc->Clone("ModelConfigNull");
134 firstPOI->setVal(poiValueForBackground);
135 mcNull->SetSnapshot(*(RooArgSet *)mcNull->GetParametersOfInterest()->snapshot());
136
137 // ----------------------------------------------------
138 // Configure a ProfileLikelihoodTestStat and a SimpleLikelihoodRatioTestStat
139 // to use simultaneously with ToyMCSampler
141 plts->SetOneSidedDiscovery(true);
142 plts->SetVarName("q_{0}/2");
143
144 // ----------------------------------------------------
145 // configure the ToyMCImportanceSampler with two test statistics
146 ToyMCSampler toymcs(*plts, 50);
147
148 // Since this tool needs to throw toy MC the PDF needs to be
149 // extended or the tool needs to know how many entries in a dataset
150 // per pseudo experiment.
151 // In the 'number counting form' where the entries in the dataset
152 // are counts, and not values of discriminating variables, the
153 // datasets typically only have one entry and the PDF is not
154 // extended.
155 if (!mc->GetPdf()->canBeExtended()) {
156 if (data->numEntries() == 1) {
157 toymcs.SetNEventsPerToy(1);
158 } else
159 cout << "Not sure what to do about this model" << endl;
160 }
161
162 // We can use PROOF to speed things along in parallel
163 // ProofConfig pc(*w, 2, "user@yourfavoriteproofcluster", false);
164 ProofConfig pc(*w, 2, "", false);
165 // toymcs.SetProofConfig(&pc); // enable proof
166
167 // instantiate the calculator
168 FrequentistCalculator freqCalc(*data, *mc, *mcNull, &toymcs);
169 freqCalc.SetToys(toys, toys); // null toys, alt toys
170
171 // Run the calculator and print result
172 HypoTestResult *freqCalcResult = freqCalc.GetHypoTest();
173 freqCalcResult->GetNullDistribution()->SetTitle("b only");
174 freqCalcResult->GetAltDistribution()->SetTitle("s+b");
175 freqCalcResult->Print();
176 double pvalue = freqCalcResult->NullPValue();
177
178 // stop timing
179 mn_t->Stop();
180 cout << "total CPU time: " << mn_t->CpuTime() << endl;
181 cout << "total real time: " << mn_t->RealTime() << endl;
182
183 // plot
184 TCanvas *c1 = new TCanvas();
185 HypoTestPlot *plot = new HypoTestPlot(*freqCalcResult, 100, -0.49, 9.51);
186 plot->SetLogYaxis(true);
187
188 // add chi2 to plot
189 int nPOI = 1;
190 TF1 *f = new TF1("f", TString::Format("1*ROOT::Math::chisquared_pdf(2*x,%d,0)", nPOI), 0, 20);
191 f->SetLineColor(kBlack);
192 f->SetLineStyle(7);
193 plot->AddTF1(f, TString::Format("#chi^{2}(2x,%d)", nPOI));
194
195 plot->Draw();
196 c1->SaveAs("standard_discovery_output.pdf");
197
198 return pvalue;
199}
#define f(i)
Definition: RSha256.hxx:104
@ kBlack
Definition: Rtypes.h:63
#define gROOT
Definition: TROOT.h:406
R__EXTERN TSystem * gSystem
Definition: TSystem.h:556
RooAbsArg * first() const
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:44
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:306
Bool_t canBeExtended() const
Definition: RooAbsPdf.h:236
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition: RooArgSet.h:126
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:261
Does a frequentist hypothesis test.
This class provides the plots for the result of a study performed with any of the HypoTestCalculatorG...
Definition: HypoTestPlot.h:22
HypoTestResult is a base class for results from hypothesis tests.
virtual Double_t NullPValue() const
Return p-value for null hypothesis.
SamplingDistribution * GetNullDistribution(void) const
void Print(const Option_t *="") const
Print out some information about the results Note: use Alt/Null labels for the hypotheses here as the...
SamplingDistribution * GetAltDistribution(void) const
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
virtual void SetSnapshot(const RooArgSet &set)
Set parameter values for a particular hypothesis if using a common PDF by saving a snapshot in the wo...
virtual ModelConfig * Clone(const char *name="") const override
clone
Definition: ModelConfig.h:54
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
Definition: ModelConfig.h:237
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
Definition: ModelConfig.h:234
ProfileLikelihoodTestStat is an implementation of the TestStatistic interface that calculates the pro...
virtual void SetVarName(const char *name)
Holds configuration options for proof and proof-lite.
Definition: ProofConfig.h:46
void AddTF1(TF1 *f, const char *title=NULL, Option_t *drawOptions="SAME")
add a TF1
void SetLogYaxis(Bool_t ly)
changes plot to log scale on y axis
void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
ToyMCSampler is an implementation of the TestStatSampler interface.
Definition: ToyMCSampler.h:72
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
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)
The Canvas class.
Definition: TCanvas.h:27
1-Dim function class
Definition: TF1.h:210
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
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:3942
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
Stopwatch class.
Definition: TStopwatch.h:28
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
Definition: TStopwatch.cxx:110
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
Definition: TStopwatch.cxx:125
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString.
Definition: TString.cxx:2311
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:1291
return c1
Definition: legend1.C:41
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Namespace for the RooStats classes.
Definition: Asimov.h:19
static constexpr double pc
Definition: file.py:1