ROOT   Reference Guide
rs401c_FeldmanCousins.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_roostats
3/// \notebook
4/// Produces an interval on the mean signal in a number counting experiment with known background using the
5/// Feldman-Cousins technique.
6///
7/// Using the RooStats FeldmanCousins tool with 200 bins
8/// it takes 1 min and the interval is [0.2625, 10.6125]
9/// with a step size of 0.075.
10/// The interval in Feldman & Cousins's original paper is [.29, 10.81] Phys.Rev.D57:3873-3889,1998.
11///
12/// \macro_image
13/// \macro_output
14/// \macro_code
15///
16/// \author Kyle Cranmer
17
18#include "RooGlobalFunc.h"
24
25#include "RooWorkspace.h"
26#include "RooDataSet.h"
27#include "RooRealVar.h"
28#include "RooConstVar.h"
30
31#include "RooDataHist.h"
32
33#include "RooPoisson.h"
34#include "RooPlot.h"
35
36#include "TCanvas.h"
37#include "TTree.h"
38#include "TH1F.h"
39#include "TMarker.h"
40#include "TStopwatch.h"
41
42#include <iostream>
43
45using namespace RooFit;
46using namespace RooStats;
47
48void rs401c_FeldmanCousins()
49{
50 // to time the macro... about 30 s
51 TStopwatch t;
52 t.Start();
53
54 // make a simple model
55 RooRealVar x("x", "", 1, 0, 50);
56 RooRealVar mu("mu", "", 2.5, 0, 15); // with a limit on mu>=0
57 RooConstVar b("b", "", 3.);
58 RooAddition mean("mean", "", RooArgList(mu, b));
59 RooPoisson pois("pois", "", x, mean);
60 RooArgSet parameters(mu);
61
62 // create a toy dataset
63 RooDataSet *data = pois.generate(RooArgSet(x), 1);
64 data->Print("v");
65
66 TCanvas *dataCanvas = new TCanvas("dataCanvas");
67 RooPlot *frame = x.frame();
68 data->plotOn(frame);
69 frame->Draw();
70 dataCanvas->Update();
71
72 RooWorkspace *w = new RooWorkspace();
73 ModelConfig modelConfig("poissonProblem", w);
74 modelConfig.SetPdf(pois);
75 modelConfig.SetParametersOfInterest(parameters);
76 modelConfig.SetObservables(RooArgSet(x));
77 w->Print();
78
79 //////// show use of Feldman-Cousins
80 RooStats::FeldmanCousins fc(*data, modelConfig);
81 fc.SetTestSize(.05); // set size of test
83 fc.FluctuateNumDataEntries(false); // number counting analysis: dataset always has 1 entry with N events observed
84 fc.SetNBins(100); // number of points to test per parameter
85
86 // use the Feldman-Cousins tool
87 PointSetInterval *interval = (PointSetInterval *)fc.GetInterval();
88
89 // make a canvas for plots
90 TCanvas *intervalCanvas = new TCanvas("intervalCanvas");
91
92 std::cout << "is this point in the interval? " << interval->IsInInterval(parameters) << std::endl;
93
94 std::cout << "interval is [" << interval->LowerLimit(mu) << ", " << interval->UpperLimit(mu) << "]" << endl;
95
96 // using 200 bins it takes 1 min and the answer is
97 // interval is [0.2625, 10.6125] with a step size of .075
98 // The interval in Feldman & Cousins's original paper is [.29, 10.81]
99 // Phys.Rev.D57:3873-3889,1998.
100
101 // No dedicated plotting class yet, so do it by hand:
102
103 RooDataHist *parameterScan = (RooDataHist *)fc.GetPointsToScan();
104 TH1F *hist = (TH1F *)parameterScan->createHistogram("mu", 30);
105 hist->Draw();
106
107 RooArgSet *tmpPoint;
108 // loop over points to test
109 for (Int_t i = 0; i < parameterScan->numEntries(); ++i) {
110 // cout << "on parameter point " << i << " out of " << parameterScan->numEntries() << endl;
111 // get a parameter point from the list of points to test.
112 tmpPoint = (RooArgSet *)parameterScan->get(i)->clone("temp");
113
114 TMarker *mark = new TMarker(tmpPoint->getRealValue("mu"), 1, 25);
115 if (interval->IsInInterval(*tmpPoint))
116 mark->SetMarkerColor(kBlue);
117 else
118 mark->SetMarkerColor(kRed);
119
120 mark->Draw("s");
121 // delete tmpPoint;
122 // delete mark;
123 }
124 t.Stop();
125 t.Print();
126}
#define b(i)
Definition: RSha256.hxx:100
int Int_t
Definition: RtypesCore.h:45
@ kRed
Definition: Rtypes.h:66
@ kBlue
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.
TH1 * createHistogram(const char *name, const RooAbsRealLValue &xvar, 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()) const
Calls createHistogram(const char *name, const RooAbsRealLValue& xvar, const RooLinkedList& argList) c...
Definition: RooAbsData.cxx:629
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooAbsData.h:193
virtual RooPlot * plotOn(RooPlot *frame, 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()) const
Definition: RooAbsData.cxx:547
RooAddition calculates the sum of a set of RooAbsReal terms, or when constructed with two sets,...
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:29
TObject * clone(const char *newname) const override
Definition: RooArgSet.h:83
RooConstVar represent a constant real-valued object.
Definition: RooConstVar.h:26
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:37
Int_t numEntries() const override
Return the number of bins.
const RooArgSet * get() const override
Get bin centre of current bin.
Definition: RooDataHist.h:74
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
A RooPlot is a plot frame and a container for graphics objects within that frame.
Definition: RooPlot.h:44
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:691
Poisson pdf.
Definition: RooPoisson.h:19
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:39
The FeldmanCousins class (like the Feldman-Cousins technique) is essentially a specific configuration...
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.
Double_t UpperLimit(RooRealVar &param)
return upper limit on a given parameter
Double_t LowerLimit(RooRealVar &param)
return lower limit on a given parameter
virtual Bool_t IsInInterval(const RooArgSet &) const
check if parameter is in the interval
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.
The Canvas class.
Definition: TCanvas.h:23
void Update() override
Definition: TCanvas.cxx:2504
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:575
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:3073
Manages Markers.
Definition: TMarker.h:22
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
Double_t x[n]
Definition: legend1.C:17
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Namespace for the RooStats classes.
Definition: Asimov.h:19
#define mark(osub)
Definition: triangle.c:1206