Logo ROOT   6.21/01
Reference Guide
rs101_limitexample.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_roostats
3 /// \notebook
4 /// Limits: number counting experiment with uncertainty on both the background rate and signal efficiency.
5 ///
6 /// The usage of a Confidence Interval Calculator to set a limit on the signal is illustrated
7 ///
8 /// \macro_image
9 /// \macro_output
10 /// \macro_code
11 ///
12 /// \author Kyle Cranmer
13 
14 #include "RooProfileLL.h"
15 #include "RooAbsPdf.h"
17 #include "RooRealVar.h"
18 #include "RooPlot.h"
19 #include "RooDataSet.h"
20 #include "RooTreeDataStore.h"
21 #include "TTree.h"
22 #include "TCanvas.h"
23 #include "TLine.h"
24 #include "TStopwatch.h"
25 
31 #include "RooStats/ConfInterval.h"
35 #include "RooStats/RooStatsUtils.h"
36 #include "RooStats/ModelConfig.h"
37 #include "RooStats/MCMCInterval.h"
41 #include "RooFitResult.h"
42 #include "TGraph2D.h"
43 
44 #include <cassert>
45 
46 // use this order for safety on library loading
47 using namespace RooFit;
48 using namespace RooStats;
49 
50 void rs101_limitexample()
51 {
52  // --------------------------------------
53  // An example of setting a limit in a number counting experiment with uncertainty on background and signal
54 
55  // to time the macro
56  TStopwatch t;
57  t.Start();
58 
59  // --------------------------------------
60  // The Model building stage
61  // --------------------------------------
62  RooWorkspace *wspace = new RooWorkspace();
63  wspace->factory("Poisson::countingModel(obs[150,0,300], "
64  "sum(s[50,0,120]*ratioSigEff[1.,0,3.],b[100]*ratioBkgEff[1.,0.,3.]))"); // counting model
65  // wspace->factory("Gaussian::sigConstraint(ratioSigEff,1,0.05)"); // 5% signal efficiency uncertainty
66  // wspace->factory("Gaussian::bkgConstraint(ratioBkgEff,1,0.1)"); // 10% background efficiency uncertainty
67  wspace->factory("Gaussian::sigConstraint(gSigEff[1,0,3],ratioSigEff,0.05)"); // 5% signal efficiency uncertainty
68  wspace->factory("Gaussian::bkgConstraint(gSigBkg[1,0,3],ratioBkgEff,0.2)"); // 10% background efficiency uncertainty
69  wspace->factory("PROD::modelWithConstraints(countingModel,sigConstraint,bkgConstraint)"); // product of terms
70  wspace->Print();
71 
72  RooAbsPdf *modelWithConstraints = wspace->pdf("modelWithConstraints"); // get the model
73  RooRealVar *obs = wspace->var("obs"); // get the observable
74  RooRealVar *s = wspace->var("s"); // get the signal we care about
75  RooRealVar *b =
76  wspace->var("b"); // get the background and set it to a constant. Uncertainty included in ratioBkgEff
77  b->setConstant();
78 
79  RooRealVar *ratioSigEff = wspace->var("ratioSigEff"); // get uncertain parameter to constrain
80  RooRealVar *ratioBkgEff = wspace->var("ratioBkgEff"); // get uncertain parameter to constrain
81  RooArgSet constrainedParams(*ratioSigEff,
82  *ratioBkgEff); // need to constrain these in the fit (should change default behavior)
83 
84  RooRealVar *gSigEff = wspace->var("gSigEff"); // global observables for signal efficiency
85  RooRealVar *gSigBkg = wspace->var("gSigBkg"); // global obs for background efficiency
86  gSigEff->setConstant();
87  gSigBkg->setConstant();
88 
89  // Create an example dataset with 160 observed events
90  obs->setVal(160.);
91  RooDataSet *data = new RooDataSet("exampleData", "exampleData", RooArgSet(*obs));
92  data->add(*obs);
93 
94  RooArgSet all(*s, *ratioBkgEff, *ratioSigEff);
95 
96  // not necessary
97  modelWithConstraints->fitTo(*data, RooFit::Constrain(RooArgSet(*ratioSigEff, *ratioBkgEff)));
98 
99  // Now let's make some confidence intervals for s, our parameter of interest
100  RooArgSet paramOfInterest(*s);
101 
102  ModelConfig modelConfig(wspace);
103  modelConfig.SetPdf(*modelWithConstraints);
104  modelConfig.SetParametersOfInterest(paramOfInterest);
105  modelConfig.SetNuisanceParameters(constrainedParams);
106  modelConfig.SetObservables(*obs);
107  modelConfig.SetGlobalObservables(RooArgSet(*gSigEff, *gSigBkg));
108  modelConfig.SetName("ModelConfig");
109  wspace->import(modelConfig);
110  wspace->import(*data);
111  wspace->SetName("w");
112  wspace->writeToFile("rs101_ws.root");
113 
114  // First, let's use a Calculator based on the Profile Likelihood Ratio
115  // ProfileLikelihoodCalculator plc(*data, *modelWithConstraints, paramOfInterest);
116  ProfileLikelihoodCalculator plc(*data, modelConfig);
117  plc.SetTestSize(.05);
118  ConfInterval *lrinterval = plc.GetInterval(); // that was easy.
119 
120  // Let's make a plot
121  TCanvas *dataCanvas = new TCanvas("dataCanvas");
122  dataCanvas->Divide(2, 1);
123 
124  dataCanvas->cd(1);
125  LikelihoodIntervalPlot plotInt((LikelihoodInterval *)lrinterval);
126  plotInt.SetTitle("Profile Likelihood Ratio and Posterior for S");
127  plotInt.Draw();
128 
129  // Second, use a Calculator based on the Feldman Cousins technique
130  FeldmanCousins fc(*data, modelConfig);
131  fc.UseAdaptiveSampling(true);
132  fc.FluctuateNumDataEntries(false); // number counting analysis: dataset always has 1 entry with N events observed
133  fc.SetNBins(100); // number of points to test per parameter
134  fc.SetTestSize(.05);
135  // fc.SaveBeltToFile(true); // optional
136  ConfInterval *fcint = NULL;
137  fcint = fc.GetInterval(); // that was easy.
138 
139  RooFitResult *fit = modelWithConstraints->fitTo(*data, Save(true));
140 
141  // Third, use a Calculator based on Markov Chain monte carlo
142  // Before configuring the calculator, let's make a ProposalFunction
143  // that will achieve a high acceptance rate
144  ProposalHelper ph;
145  ph.SetVariables((RooArgSet &)fit->floatParsFinal());
146  ph.SetCovMatrix(fit->covarianceMatrix());
148  ph.SetCacheSize(100);
149  ProposalFunction *pdfProp = ph.GetProposalFunction(); // that was easy
150 
151  MCMCCalculator mc(*data, modelConfig);
152  mc.SetNumIters(20000); // steps to propose in the chain
153  mc.SetTestSize(.05); // 95% CL
154  mc.SetNumBurnInSteps(40); // ignore first N steps in chain as "burn in"
155  mc.SetProposalFunction(*pdfProp);
156  mc.SetLeftSideTailFraction(0.5); // find a "central" interval
157  MCMCInterval *mcInt = (MCMCInterval *)mc.GetInterval(); // that was easy
158 
159  // Get Lower and Upper limits from Profile Calculator
160  cout << "Profile lower limit on s = " << ((LikelihoodInterval *)lrinterval)->LowerLimit(*s) << endl;
161  cout << "Profile upper limit on s = " << ((LikelihoodInterval *)lrinterval)->UpperLimit(*s) << endl;
162 
163  // Get Lower and Upper limits from FeldmanCousins with profile construction
164  if (fcint != NULL) {
165  double fcul = ((PointSetInterval *)fcint)->UpperLimit(*s);
166  double fcll = ((PointSetInterval *)fcint)->LowerLimit(*s);
167  cout << "FC lower limit on s = " << fcll << endl;
168  cout << "FC upper limit on s = " << fcul << endl;
169  TLine *fcllLine = new TLine(fcll, 0, fcll, 1);
170  TLine *fculLine = new TLine(fcul, 0, fcul, 1);
171  fcllLine->SetLineColor(kRed);
172  fculLine->SetLineColor(kRed);
173  fcllLine->Draw("same");
174  fculLine->Draw("same");
175  dataCanvas->Update();
176  }
177 
178  // Plot MCMC interval and print some statistics
179  MCMCIntervalPlot mcPlot(*mcInt);
180  mcPlot.SetLineColor(kMagenta);
181  mcPlot.SetLineWidth(2);
182  mcPlot.Draw("same");
183 
184  double mcul = mcInt->UpperLimit(*s);
185  double mcll = mcInt->LowerLimit(*s);
186  cout << "MCMC lower limit on s = " << mcll << endl;
187  cout << "MCMC upper limit on s = " << mcul << endl;
188  cout << "MCMC Actual confidence level: " << mcInt->GetActualConfidenceLevel() << endl;
189 
190  // 3-d plot of the parameter points
191  dataCanvas->cd(2);
192  // also plot the points in the markov chain
193  RooDataSet *chainData = mcInt->GetChainAsDataSet();
194 
195  assert(chainData);
196  std::cout << "plotting the chain data - nentries = " << chainData->numEntries() << std::endl;
197  TTree *chain = RooStats::GetAsTTree("chainTreeData", "chainTreeData", *chainData);
198  assert(chain);
199  chain->SetMarkerStyle(6);
200  chain->SetMarkerColor(kRed);
201 
202  chain->Draw("s:ratioSigEff:ratioBkgEff", "nll_MarkovChain_local_", "box"); // 3-d box proportional to posterior
203 
204  // the points used in the profile construction
205  RooDataSet *parScanData = (RooDataSet *)fc.GetPointsToScan();
206  assert(parScanData);
207  std::cout << "plotting the scanned points used in the frequentist construction - npoints = "
208  << parScanData->numEntries() << std::endl;
209  // getting the tree and drawing it -crashes (very strange....);
210  // TTree* parameterScan = RooStats::GetAsTTree("parScanTreeData","parScanTreeData",*parScanData);
211  // assert(parameterScan);
212  // parameterScan->Draw("s:ratioSigEff:ratioBkgEff","","goff");
213  TGraph2D *gr = new TGraph2D(parScanData->numEntries());
214  for (int ievt = 0; ievt < parScanData->numEntries(); ++ievt) {
215  const RooArgSet *evt = parScanData->get(ievt);
216  double x = evt->getRealValue("ratioBkgEff");
217  double y = evt->getRealValue("ratioSigEff");
218  double z = evt->getRealValue("s");
219  gr->SetPoint(ievt, x, y, z);
220  // std::cout << ievt << " " << x << " " << y << " " << z << std::endl;
221  }
222  gr->SetMarkerStyle(24);
223  gr->Draw("P SAME");
224 
225  delete wspace;
226  delete lrinterval;
227  delete mcInt;
228  delete fcint;
229  delete data;
230 
231  // print timing info
232  t.Stop();
233  t.Print();
234 }
virtual const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event &#39;index&#39;.
ProposalFunction is an interface for all proposal functions that would be used with a Markov Chain Mo...
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
const RooArgList & floatParsFinal() const
Definition: RooFitResult.h:110
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
This class provides simple and straightforward utilities to plot a MCMCInterval object.
void Print(Option_t *option="") const
Print the real and cpu time passed between the start and stop events.
Definition: TStopwatch.cxx:219
LikelihoodInterval is a concrete implementation of the RooStats::ConfInterval interface.
Definition: Rtypes.h:64
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
const TMatrixDSym & covarianceMatrix() const
Return covariance matrix.
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition: TNamed.cxx:140
virtual ProposalFunction * GetProposalFunction()
virtual void SetUpdateProposalParameters(Bool_t updateParams)
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:696
This class provides simple and straightforward utilities to plot a LikelihoodInterval object...
The ProfileLikelihoodCalculator is a concrete implementation of CombinedCalculator (the interface cla...
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
Definition: RooFitResult.h:40
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition: TObject.cxx:195
TTree * GetAsTTree(TString name, TString desc, const RooDataSet &data)
virtual void Draw(Option_t *chopt="")
Draw this graph with its current attributes.
Definition: TGraph.cxx:753
static struct mg_connection * fc(struct mg_context *ctx)
Definition: civetweb.c:3728
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
Double_t x[n]
Definition: legend1.C:17
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
Definition: TAttMarker.h:38
static constexpr double s
RooRealVar represents a fundamental (non-derived) real-valued object.
Definition: RooRealVar.h:36
virtual void setVal(Double_t value)
Set value of variable to &#39;value&#39;.
Definition: RooRealVar.cxx:252
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
void setConstant(Bool_t value=kTRUE)
virtual void SetCacheSize(Int_t size)
A simple line.
Definition: TLine.h:23
Bool_t writeToFile(const char *fileName, Bool_t recreate=kTRUE)
Save this current workspace into given file.
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition: TAttMarker.h:40
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:31
TGraphErrors * gr
Definition: legend1.C:25
The FeldmanCousins class (like the Feldman-Cousins technique) is essentially a specific configuration...
The Canvas class.
Definition: TCanvas.h:31
ConfInterval is an interface class for a generic interval in the RooStats framework.
Definition: ConfInterval.h:35
Namespace for the RooStats classes.
Definition: Asimov.h:20
PointSetInterval is a concrete implementation of the ConfInterval interface.
RooAbsPdf * pdf(const char *name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
virtual Double_t LowerLimit(RooRealVar &param)
get the lowest value of param that is within the confidence interval
virtual void Draw(Option_t *opt)
Default Draw method for all objects.
Definition: TTree.h:411
Double_t y[n]
Definition: legend1.C:17
virtual void SetCovMatrix(const TMatrixDSym &covMatrix)
RooRealVar * var(const char *name) const
Retrieve real-valued variable (RooRealVar) with given name. A null pointer is returned if not found...
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
RooFactoryWSTool & factory()
Return instance to factory tool.
RooCmdArg Save(Bool_t flag=kTRUE)
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition: TGraph.cxx:2257
RooAbsPdf, the base class of all PDFs
Definition: RooAbsPdf.h:40
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)
Automatic pad generation by division.
Definition: TPad.cxx:1165
Bool_t import(const RooAbsArg &arg, const RooCmdArg &arg1=RooCmdArg(), const RooCmdArg &arg2=RooCmdArg(), const RooCmdArg &arg3=RooCmdArg(), const RooCmdArg &arg4=RooCmdArg(), const RooCmdArg &arg5=RooCmdArg(), const RooCmdArg &arg6=RooCmdArg(), const RooCmdArg &arg7=RooCmdArg(), const RooCmdArg &arg8=RooCmdArg(), const RooCmdArg &arg9=RooCmdArg())
Import a RooAbsArg object, e.g.
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:472
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
virtual RooFitResult * fitTo(RooAbsData &data, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none())
Fit PDF to given dataset.
Definition: RooAbsPdf.cxx:1256
A TTree represents a columnar dataset.
Definition: TTree.h:72
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
Definition: MCMCInterval.h:30
virtual void SetVariables(RooArgList &vars)
virtual void Update()
Update canvas pad buffers.
Definition: TCanvas.cxx:2339
Graphics object made of three arrays X, Y and Z with the same number of points each.
Definition: TGraph2D.h:40
Bayesian Calculator estimating an interval or a credible region using the Markov-Chain Monte Carlo me...
void Print(Option_t *opts=0) const
Print contents of the workspace.
RooCmdArg Constrain(const RooArgSet &params)
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:306
virtual void add(const RooArgSet &row, Double_t weight=1.0, Double_t weightError=0) override
Add a data point, with its coordinates specified in the &#39;data&#39; argset, to the data set...
Stopwatch class.
Definition: TStopwatch.h:28