Logo ROOT   6.21/01
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 
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 }
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:1279
Holds configuration options for proof and proof-lite.
Definition: ProofConfig.h:46
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
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
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:3925
void Print(const Option_t *="") const
Print out some information about the results Note: use Alt/Null labels for the hypotheses here as the...
return c1
Definition: legend1.C:41
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
Definition: Rtypes.h:63
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format...
Definition: TFile.h:48
HypoTestResult is a base class for results from hypothesis tests.
#define gROOT
Definition: TROOT.h:415
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
#define f(i)
Definition: RSha256.hxx:104
This class provides the plots for the result of a study performed with any of the HypoTestCalculatorG...
Definition: HypoTestPlot.h:22
SamplingDistribution * GetAltDistribution(void) const
virtual ModelConfig * Clone(const char *name="") const override
clone
Definition: ModelConfig.h:54
void SetLogYaxis(Bool_t ly)
changes plot to log scale on y axis
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
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
RooRealVar represents a fundamental (non-derived) real-valued object.
Definition: RooRealVar.h:36
RooAbsData * data(const char *name) const
Retrieve dataset (binned or unbinned) with given name. A null pointer is returned if not found...
virtual void setVal(Double_t value)
Set value of variable to &#39;value&#39;.
Definition: RooRealVar.cxx:252
R__EXTERN TSystem * gSystem
Definition: TSystem.h:557
static constexpr double pc
RooAbsArg * first() const
virtual Double_t NullPValue() const
Return p-value for null hypothesis.
SamplingDistribution * GetNullDistribution(void) const
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:44
ProfileLikelihoodTestStat is an implementation of the TestStatistic interface that calculates the pro...
ToyMCSampler is an implementation of the TestStatSampler interface.
Definition: ToyMCSampler.h:71
Bool_t canBeExtended() const
Definition: RooAbsPdf.h:236
TObject * obj(const char *name) const
Return any type of object (RooAbsArg, RooAbsData or generic object) with given name) ...
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition: RooArgSet.h:134
The Canvas class.
Definition: TCanvas.h:31
Namespace for the RooStats classes.
Definition: Asimov.h:20
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
Definition: ModelConfig.h:234
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
Definition: ModelConfig.h:237
void AddTF1(TF1 *f, const char *title=NULL, Option_t *drawOptions="SAME")
add a TF1
Definition: file.py:1
1-Dim function class
Definition: TF1.h:211
Does a frequentist hypothesis test.
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
void Print(Option_t *opts=0) const
Print contents of the workspace.
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 void SetVarName(const char *name)
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:306
Stopwatch class.
Definition: TStopwatch.h:28