Logo ROOT   6.08/07
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 
56 double StandardFrequentistDiscovery(
57  const char* infile = "",
58  const char* workspaceName = "channel1",
59  const char* modelConfigNameSB = "ModelConfig",
60  const char* dataName = "obsData",
61  int toys = 1000,
62  double poiValueForBackground = 0.0,
63  double poiValueForSignal = 1.0
64 ) {
65 
66  // The workspace contains the model for s+b. The b model is "autogenerated"
67  // by copying s+b and setting the one parameter of interest to zero.
68  // To keep the script simple, multiple parameters of interest or different
69  // functional forms of the b model are not supported.
70 
71  // for now, assume there is only one parameter of interest, and these are
72  // its values:
73 
74  // -------------------------------------------------------
75  // First part is just to access a user-defined file
76  // or create the standard example file if it doesn't exist
77  const char* filename = "";
78  if (!strcmp(infile,"")) {
79  filename = "results/example_channel1_GammaExample_model.root";
80  bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code
81  // if file does not exists generate with histfactory
82  if (!fileExist) {
83 #ifdef _WIN32
84  cout << "HistFactory file cannot be generated on Windows - exit" << endl;
85  return -1;
86 #endif
87  // Normally this would be run on the command line
88  cout <<"will run standard hist2workspace example"<<endl;
89  gROOT->ProcessLine(".! prepareHistFactory .");
90  gROOT->ProcessLine(".! hist2workspace config/example.xml");
91  cout <<"\n\n---------------------"<<endl;
92  cout <<"Done creating example input"<<endl;
93  cout <<"---------------------\n\n"<<endl;
94  }
95 
96  }
97  else
98  filename = infile;
99 
100  // Try to open the file
101  TFile *file = TFile::Open(filename);
102 
103  // if input file was specified byt not found, quit
104  if(!file ){
105  cout <<"StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
106  return -1;
107  }
108 
109 
110  // -------------------------------------------------------
111  // Tutorial starts here
112  // -------------------------------------------------------
113 
114  TStopwatch *mn_t = new TStopwatch;
115  mn_t->Start();
116 
117  // get the workspace out of the file
118  RooWorkspace* w = (RooWorkspace*) file->Get(workspaceName);
119  if (!w) {
120  cout << "workspace not found" << endl;
121  return -1.0;
122  }
123 
124  // get the modelConfig out of the file
125  ModelConfig* mc = (ModelConfig*) w->obj(modelConfigNameSB);
126 
127  // get the data out of the file
128  RooAbsData* data = w->data(dataName);
129 
130  // make sure ingredients are found
131  if (!data || !mc) {
132  w->Print();
133  cout << "data or ModelConfig was not found" << endl;
134  return -1.0;
135  }
136 
137 
138  RooRealVar* firstPOI = (RooRealVar*) mc->GetParametersOfInterest()->first();
139  firstPOI->setVal(poiValueForSignal);
141  // create null model
142  ModelConfig *mcNull = mc->Clone("ModelConfigNull");
143  firstPOI->setVal(poiValueForBackground);
144  mcNull->SetSnapshot(*(RooArgSet*)mcNull->GetParametersOfInterest()->snapshot());
145 
146 
147 
148  // ----------------------------------------------------
149  // Configure a ProfileLikelihoodTestStat and a SimpleLikelihoodRatioTestStat
150  // to use simultaneously with ToyMCSampler
152  plts->SetOneSidedDiscovery(true);
153  plts->SetVarName( "q_{0}/2" );
154 
155  // ----------------------------------------------------
156  // configure the ToyMCImportanceSampler with two test statistics
157  ToyMCSampler toymcs(*plts, 50);
158 
159 
160 
161  // Since this tool needs to throw toy MC the PDF needs to be
162  // extended or the tool needs to know how many entries in a dataset
163  // per pseudo experiment.
164  // In the 'number counting form' where the entries in the dataset
165  // are counts, and not values of discriminating variables, the
166  // datasets typically only have one entry and the PDF is not
167  // extended.
168  if (!mc->GetPdf()->canBeExtended()) {
169  if (data->numEntries() == 1) {
170  toymcs.SetNEventsPerToy(1);
171  } else cout << "Not sure what to do about this model" << endl;
172  }
173 
174  // We can use PROOF to speed things along in parallel
175  // ProofConfig pc(*w, 2, "user@yourfavoriteproofcluster", false);
176  ProofConfig pc(*w, 2, "", false);
177  //toymcs.SetProofConfig(&pc); // enable proof
178 
179 
180  // instantiate the calculator
181  FrequentistCalculator freqCalc(*data, *mc, *mcNull, &toymcs);
182  freqCalc.SetToys( toys,toys ); // null toys, alt toys
183 
184  // Run the calculator and print result
185  HypoTestResult* freqCalcResult = freqCalc.GetHypoTest();
186  freqCalcResult->GetNullDistribution()->SetTitle( "b only" );
187  freqCalcResult->GetAltDistribution()->SetTitle( "s+b" );
188  freqCalcResult->Print();
189  double pvalue = freqCalcResult->NullPValue();
190 
191  // stop timing
192  mn_t->Stop();
193  cout << "total CPU time: " << mn_t->CpuTime() << endl;
194  cout << "total real time: " << mn_t->RealTime() << endl;
195 
196  // plot
197  TCanvas* c1 = new TCanvas();
198  HypoTestPlot *plot = new HypoTestPlot(*freqCalcResult, 100, -0.49, 9.51 );
199  plot->SetLogYaxis(true);
200 
201  // add chi2 to plot
202  int nPOI = 1;
203  TF1* f = new TF1("f", TString::Format("1*ROOT::Math::chisquared_pdf(2*x,%d,0)",nPOI), 0,20);
204  f->SetLineColor( kBlack );
205  f->SetLineStyle( 7 );
206  plot->AddTF1( f, TString::Format("#chi^{2}(2x,%d)",nPOI) );
207 
208  plot->Draw();
209  c1->SaveAs("standard_discovery_output.pdf");
210 
211 
212  return pvalue;
213 }
214 
215 
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:1266
Holds configuration options for proof and proof-lite.
Definition: ProofConfig.h:49
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:52
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
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
Definition: Rtypes.h:60
HypoTestResult is a base class for results from hypothesis tests.
#define gROOT
Definition: TROOT.h:364
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
This class provides the plots for the result of a study performed with any of the HypoTestCalculatorG...
Definition: HypoTestPlot.h:37
SamplingDistribution * GetAltDistribution(void) const
void SetLogYaxis(Bool_t ly)
changes plot to log scale on y axis
virtual ModelConfig * Clone(const char *name="") const
clone
Definition: ModelConfig.h:76
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=1, Int_t netopt=0)
Create / open a file.
Definition: TFile.cxx:3907
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:2335
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:37
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:205
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:46
RooAbsCollection * snapshot(Bool_t deepCopy=kTRUE) const
Take a snap shot of current collection contents: An owning collection is returned containing clones o...
R__EXTERN TSystem * gSystem
Definition: TSystem.h:549
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:37
ProfileLikelihoodTestStat is an implementation of the TestStatistic interface that calculates the pro...
ToyMCSampler is an implementation of the TestStatSampler interface.
Definition: ToyMCSampler.h:99
Bool_t canBeExtended() const
Definition: RooAbsPdf.h:216
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:41
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:244
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
Definition: ModelConfig.h:247
void AddTF1(TF1 *f, const char *title=NULL, Option_t *drawOptions="SAME")
add a TF1
double f(double x)
virtual void SetLineStyle(Style_t lstyle)
Set the line style.
Definition: TAttLine.h:48
Definition: file.py:1
1-Dim function class
Definition: TF1.h:149
Hypothesis Test Calculator using a full frequentist procedure for sampling the test statistic distrib...
virtual void SaveAs(const char *filename="", Option_t *option="") const
Save Pad contents in a file in one of various formats.
Definition: TPad.cxx:5037
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:155
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:42
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:269
Stopwatch class.
Definition: TStopwatch.h:30