Logo ROOT  
Reference Guide
StandardFrequentistDiscovery.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_roostats
3 /// \notebook
4 /// \brief 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 
36 #include "RooStats/ModelConfig.h"
39 #include "RooStats/HypoTestPlot.h"
46 
48 #include "TSystem.h"
49 
50 #include <vector>
51 
52 using namespace RooFit;
53 using namespace RooStats;
54 
55 double 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 }
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
RooRealVar::setVal
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:216
TH1F.h
f
#define f(i)
Definition: RSha256.hxx:104
RooAbsCollection::first
RooAbsArg * first() const
Definition: RooAbsCollection.h:176
RooAbsData
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:46
RooNumIntConfig.h
RooStats::HypoTestResult::GetNullDistribution
SamplingDistribution * GetNullDistribution(void) const
Definition: HypoTestResult.h:70
RooStats::HypoTestResult::Print
void Print(const Option_t *="") const
Print out some information about the results Note: use Alt/Null labels for the hypotheses here as the...
Definition: HypoTestResult.cxx:339
RooStats::HypoTestResult::NullPValue
virtual Double_t NullPValue() const
Return p-value for null hypothesis.
Definition: HypoTestResult.h:45
TStopwatch.h
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:3995
TGeant4Unit::pc
static constexpr double pc
Definition: TGeant4SystemOfUnits.h:130
RooStats::ModelConfig::Clone
virtual ModelConfig * Clone(const char *name="") const override
clone
Definition: ModelConfig.h:54
TString::Format
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
TCanvas.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
HypoTestResult.h
ToyMCImportanceSampler.h
TFile.h
RooStats::HypoTestPlot
This class provides the plots for the result of a study performed with any of the HypoTestCalculatorG...
Definition: HypoTestPlot.h:22
TROOT.h
TStopwatch::RealTime
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
RooArgSet::snapshot
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition: RooArgSet.h:134
kBlack
@ kBlack
Definition: Rtypes.h:65
RooStats::FrequentistCalculator
Does a frequentist hypothesis test.
Definition: FrequentistCalculator.h:25
LikelihoodInterval.h
TSystem.h
FrequentistCalculator.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
SimpleLikelihoodRatioTestStat.h
RooStats::ProofConfig
Holds configuration options for proof and proof-lite.
Definition: ProofConfig.h:46
RooRandom.h
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
TStopwatch::Start
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
TNamed::SetTitle
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
TStopwatch::CpuTime
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
RooStats::HypoTestResult::GetAltDistribution
SamplingDistribution * GetAltDistribution(void) const
Definition: HypoTestResult.h:71
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
SamplingDistribution.h
gSystem
R__EXTERN TSystem * gSystem
Definition: TSystem.h:559
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
TF1.h
RooStats::SamplingDistPlot::Draw
void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: SamplingDistPlot.cxx:297
TCanvas
The Canvas class.
Definition: TCanvas.h:23
RooStats::ProfileLikelihoodTestStat::SetOneSidedDiscovery
void SetOneSidedDiscovery(Bool_t flag=true)
Definition: ProfileLikelihoodTestStat.h:88
RooStats
Namespace for the RooStats classes.
Definition: Asimov.h:19
file
Definition: file.py:1
RooStats::ToyMCSampler
ToyMCSampler is an implementation of the TestStatSampler interface.
Definition: ToyMCSampler.h:73
TStopwatch
Stopwatch class.
Definition: TStopwatch.h:28
RooStats::ProfileLikelihoodTestStat::SetVarName
virtual void SetVarName(const char *name)
Definition: ProfileLikelihoodTestStat.h:137
ProfileLikelihoodCalculator.h
RooStats::HypoTestResult
HypoTestResult is a base class for results from hypothesis tests.
Definition: HypoTestResult.h:22
RooStats::ProfileLikelihoodTestStat
ProfileLikelihoodTestStat is an implementation of the TestStatistic interface that calculates the pro...
Definition: ProfileLikelihoodTestStat.h:32
LikelihoodIntervalPlot.h
TF1
1-Dim function class
Definition: TF1.h:213
TStopwatch::Stop
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
RooRealSumPdf.h
RooRealVar
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:36
RooStats::ModelConfig
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
HypoTestPlot.h
RooStats::ModelConfig::SetSnapshot
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...
Definition: ModelConfig.cxx:209
RooAbsPdf::canBeExtended
Bool_t canBeExtended() const
Definition: RooAbsPdf.h:236
RooStats::SamplingDistPlot::AddTF1
void AddTF1(TF1 *f, const char *title=NULL, Option_t *drawOptions="SAME")
add a TF1
Definition: SamplingDistPlot.cxx:235
RooArgSet
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
gROOT
#define gROOT
Definition: TROOT.h:406
ProfileLikelihoodTestStat.h
c1
return c1
Definition: legend1.C:41
RooStats::SamplingDistPlot::SetLogYaxis
void SetLogYaxis(Bool_t ly)
changes plot to log scale on y axis
Definition: SamplingDistPlot.h:93