ROOT   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
41#include "RooFitResult.h"
42#include "TGraph2D.h"
43
44#include <cassert>
45
47using namespace RooFit;
48using namespace RooStats;
49
50void 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));
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);
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
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}
#define b(i)
Definition: RSha256.hxx:100
@ kRed
Definition: Rtypes.h:66
@ kMagenta
Definition: Rtypes.h:66
static struct mg_connection * fc(struct mg_context *ctx)
Definition: civetweb.c:3728
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.
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Definition: RooAbsData.cxx:307
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:1410
void setConstant(Bool_t value=kTRUE)
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:29
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
virtual const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
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 'data' argset, to the data set.
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
Definition: RooFitResult.h:40
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:39
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:282
ConfInterval is an interface class for a generic interval in the RooStats framework.
Definition: ConfInterval.h:35
The FeldmanCousins class (like the Feldman-Cousins technique) is essentially a specific configuration...
This class provides simple and straightforward utilities to plot a LikelihoodInterval object.
LikelihoodInterval is a concrete implementation of the RooStats::ConfInterval interface.
Bayesian Calculator estimating an interval or a credible region using the Markov-Chain Monte Carlo me...
This class provides simple and straightforward utilities to plot a MCMCInterval object.
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
Definition: MCMCInterval.h:33
virtual Double_t LowerLimit(RooRealVar &param)
get the lowest value of param that is within the confidence interval
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
PointSetInterval is a concrete implementation of the ConfInterval interface.
The ProfileLikelihoodCalculator is a concrete implementation of CombinedCalculator (the interface cla...
ProposalFunction is an interface for all proposal functions that would be used with a Markov Chain Mo...
virtual void SetCovMatrix(const TMatrixDSym &covMatrix)
virtual void SetUpdateProposalParameters(Bool_t updateParams)
virtual ProposalFunction * GetProposalFunction()
virtual void SetVariables(RooArgList &vars)
virtual void SetCacheSize(Int_t size)
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
void Print(Option_t *opts=0) const
Print contents of the workspace.
Bool_t writeToFile(const char *fileName, Bool_t recreate=kTRUE)
Save this current workspace into given file.
RooRealVar * var(const char *name) const
Retrieve real-valued variable (RooRealVar) with given name. A null pointer is returned if not found.
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.
RooFactoryWSTool & factory()
Return instance to factory tool.
RooAbsPdf * pdf(const char *name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
Definition: TAttMarker.h:38
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition: TAttMarker.h:40
The Canvas class.
Definition: TCanvas.h:23
Definition: TCanvas.cxx:708
void Update() override
Definition: TCanvas.cxx:2504
Graphics object made of three arrays X, Y and Z with the same number of points each.
Definition: TGraph2D.h:41
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition: TGraph.cxx:2284
virtual void Draw(Option_t *chopt="")
Draw this graph with its current attributes.
Definition: TGraph.cxx:769
A simple line.
Definition: TLine.h:22
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition: TNamed.cxx:140
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition: TObject.cxx:197
void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0) override
Stopwatch class.
Definition: TStopwatch.h:28
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
void Print(Option_t *option="") const
Print the real and cpu time passed between the start and stop events.
Definition: TStopwatch.cxx:219
A TTree represents a columnar dataset.
Definition: TTree.h:79
virtual void Draw(Option_t *opt)
Default Draw method for all objects.
Definition: TTree.h:428
RooCmdArg Constrain(const RooArgSet &params)
RooCmdArg Save(Bool_t flag=kTRUE)
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
TGraphErrors * gr
Definition: legend1.C:25