ROOT   Reference Guide
rs401c_FeldmanCousins.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_roostats
3 /// \notebook
4 /// \brief 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"
19 #include "RooStats/ConfInterval.h"
23 #include "RooStats/ModelConfig.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
45 using namespace RooFit;
46 using namespace RooStats;
47
48 void 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 }
mark
#define mark(osub)
Definition: triangle.c:1206
RooWorkspace.h
RooPlot::Draw
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:691
TH1F.h
PointSetInterval.h
RooStats::PointSetInterval::UpperLimit
Double_t UpperLimit(RooRealVar &param)
return upper limit on a given parameter
Definition: PointSetInterval.cxx:147
RooArgSet::getRealValue
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:474
fc
static struct mg_connection * fc(struct mg_context *ctx)
Definition: civetweb.c:3728
ConfidenceBelt.h
TStopwatch.h
RooArgList
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
RooStats::FeldmanCousins
The FeldmanCousins class (like the Feldman-Cousins technique) is essentially a specific configuration...
Definition: FeldmanCousins.h:33
RooDataHist::get
virtual const RooArgSet * get() const
Definition: RooDataHist.h:78
RooAbsData::Print
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooAbsData.h:182
ConfInterval.h
x
Double_t x[n]
Definition: legend1.C:17
RooConstVar
RooConstVar represent a constant real-valued object.
Definition: RooConstVar.h:25
TCanvas.h
TTree.h
RooDataSet.h
b
#define b(i)
Definition: RSha256.hxx:100
TStopwatch::Print
void Print(Option_t *option="") const
Print the real and cpu time passed between the start and stop events.
Definition: TStopwatch.cxx:219
RooDataHist
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:39
RooWorkspace::Print
void Print(Option_t *opts=0) const
Print contents of the workspace.
Definition: RooWorkspace.cxx:2194
TMarker.h
ModelConfig.h
RooFit
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition: RooCFunction1Binding.h:29
RooAbsData::plotOn
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
RooStats::PointSetInterval::IsInInterval
virtual Bool_t IsInInterval(const RooArgSet &) const
check if parameter is in the interval
Definition: PointSetInterval.cxx:76
RooDataHist.h
RooAbsData::createHistogram
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
RooPlot.h
RooAddition calculates the sum of a set of RooAbsReal terms, or when constructed with two sets,...
RooPlot
A RooPlot is a plot frame and a container for graphics objects within that frame.
Definition: RooPlot.h:44
RooPoisson
Poisson pdf.
Definition: RooPoisson.h:19
TStopwatch::Start
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
RooRealVar.h
kRed
@ kRed
Definition: Rtypes.h:66
RooConstVar.h
RooGlobalFunc.h
FeldmanCousins.h
RooArgSet::clone
virtual TObject * clone(const char *newname) const
Definition: RooArgSet.h:84
RooPoisson.h
RooWorkspace
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
TCanvas
The Canvas class.
Definition: TCanvas.h:23
RooDataHist::numEntries
virtual Int_t numEntries() const
Return the number of bins.
Definition: RooDataHist.cxx:1694
RooStats
Namespace for the RooStats classes.
Definition: Asimov.h:19
TH1F
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:573
TStopwatch
Stopwatch class.
Definition: TStopwatch.h:28
TCanvas::Update
virtual void Update()
Definition: TCanvas.cxx:2500
kBlue
@ kBlue
Definition: Rtypes.h:66
RooDataSet
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
TMarker
Manages Markers.
Definition: TMarker.h:22
TStopwatch::Stop
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
RooRealVar
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:36
RooStats::PointSetInterval
PointSetInterval is a concrete implementation of the ConfInterval interface.
Definition: PointSetInterval.h:21
RooStats::ModelConfig
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
RooStats::PointSetInterval::LowerLimit
Double_t LowerLimit(RooRealVar &param)
return lower limit on a given parameter
Definition: PointSetInterval.cxx:160
RooArgSet
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
int
TH1::Draw
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:3050