Logo ROOT   6.21/01
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"
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"
29 #include "RooAddition.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 
44 // use this order for safety on library loading
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
82  fc.UseAdaptiveSampling(true);
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 }
Poisson pdf.
Definition: RooPoisson.h:19
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
virtual const RooArgSet * get() const
Definition: RooDataHist.h:79
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
Calls RooPlot* plotOn(RooPlot* frame, const RooLinkedList& cmdList) const ;.
Definition: RooAbsData.cxx:549
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
void Print(Option_t *option="") const
Print the real and cpu time passed between the start and stop events.
Definition: TStopwatch.cxx:219
Definition: Rtypes.h:64
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
#define mark(osub)
Definition: triangle.c:1206
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:571
Manages Markers.
Definition: TMarker.h:23
int Int_t
Definition: RtypesCore.h:41
static struct mg_connection * fc(struct mg_context *ctx)
Definition: civetweb.c:3728
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:40
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...
RooConstVar represent a constant real-valued object.
Definition: RooConstVar.h:25
virtual TObject * clone(const char *newname) const
Definition: RooArgSet.h:84
RooRealVar represents a fundamental (non-derived) real-valued object.
Definition: RooRealVar.h:36
Double_t LowerLimit(RooRealVar &param)
return lower limit on a given parameter
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2981
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:31
Double_t UpperLimit(RooRealVar &param)
return upper limit on a given parameter
A RooPlot is a plot frame and a container for graphics objects within that frame. ...
Definition: RooPlot.h:44
The FeldmanCousins class (like the Feldman-Cousins technique) is essentially a specific configuration...
The Canvas class.
Definition: TCanvas.h:31
Namespace for the RooStats classes.
Definition: Asimov.h:20
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooAbsData.h:171
PointSetInterval is a concrete implementation of the ConfInterval interface.
virtual Int_t numEntries() const
Return the number of bins.
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
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
Definition: Rtypes.h:64
RooAddition calculates the sum of a set of RooAbsReal terms, or when constructed with two sets...
Definition: RooAddition.h:26
virtual void Update()
Update canvas pad buffers.
Definition: TCanvas.cxx:2339
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:631
void Print(Option_t *opts=0) const
Print contents of the workspace.
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
virtual Bool_t IsInInterval(const RooArgSet &) const
check if parameter is in the interval
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:712
Stopwatch class.
Definition: TStopwatch.h:28