# Rs 4 0 1C_ Feldman Cousins¶

Produces an interval on the mean signal in a number counting experiment with known background using the Feldman-Cousins technique.

Using the RooStats FeldmanCousins tool with 200 bins it takes 1 min and the interval is [0.2625, 10.6125] with a step size of 0.075. The interval in Feldman & Cousins's original paper is [.29, 10.81] Phys.Rev.D57:3873-3889,1998.

Author: Kyle Cranmer
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Thursday, June 30, 2022 at 09:51 AM.

In [ ]:
%%cpp -d
#include "RooGlobalFunc.h"
#include "RooStats/ConfInterval.h"
#include "RooStats/PointSetInterval.h"
#include "RooStats/ConfidenceBelt.h"
#include "RooStats/FeldmanCousins.h"
#include "RooStats/ModelConfig.h"

#include "RooWorkspace.h"
#include "RooDataSet.h"
#include "RooRealVar.h"
#include "RooConstVar.h"

#include "RooDataHist.h"

#include "RooPoisson.h"
#include "RooPlot.h"

#include "TCanvas.h"
#include "TTree.h"
#include "TH1F.h"
#include "TMarker.h"
#include "TStopwatch.h"

#include <iostream>


In [ ]:
%%cpp -d
// This is a workaround to make sure the namespace is used inside functions
using namespace RooFit;
using namespace RooStats;


To time the macro... about 30 s

In [ ]:
TStopwatch t;
t.Start();


Make a simple model

In [ ]:
RooRealVar x("x", "", 1, 0, 50);
RooRealVar mu("mu", "", 2.5, 0, 15); // with a limit on mu>=0
RooConstVar b("b", "", 3.);
RooPoisson pois("pois", "", x, mean);
RooArgSet parameters(mu);


Create a toy dataset

In [ ]:
RooDataSet *data = pois.generate(RooArgSet(x), 1);
data->Print("v");

TCanvas *dataCanvas = new TCanvas("dataCanvas");
RooPlot *frame = x.frame();
data->plotOn(frame);
frame->Draw();
dataCanvas->Update();

RooWorkspace *w = new RooWorkspace();
ModelConfig modelConfig("poissonProblem", w);
modelConfig.SetPdf(pois);
modelConfig.SetParametersOfInterest(parameters);
modelConfig.SetObservables(RooArgSet(x));
w->Print();


////// show use of feldman-cousins

In [ ]:
RooStats::FeldmanCousins fc(*data, modelConfig);
fc.SetTestSize(.05); // set size of test
fc.FluctuateNumDataEntries(false); // number counting analysis: dataset always has 1 entry with N events observed
fc.SetNBins(100);                  // number of points to test per parameter


Use the feldman-cousins tool

In [ ]:
PointSetInterval *interval = (PointSetInterval *)fc.GetInterval();


Make a canvas for plots

In [ ]:
TCanvas *intervalCanvas = new TCanvas("intervalCanvas");

std::cout << "is this point in the interval? " << interval->IsInInterval(parameters) << std::endl;

std::cout << "interval is [" << interval->LowerLimit(mu) << ", " << interval->UpperLimit(mu) << "]" << endl;


Using 200 bins it takes 1 min and the answer is interval is [0.2625, 10.6125] with a step size of .075 The interval in Feldman & Cousins's original paper is [.29, 10.81] Phys.Rev.D57:3873-3889,1998.

No dedicated plotting class yet, so do it by hand:

In [ ]:
RooDataHist *parameterScan = (RooDataHist *)fc.GetPointsToScan();
TH1F *hist = (TH1F *)parameterScan->createHistogram("mu", Binning(30));
hist->Draw();

RooArgSet *tmpPoint;


Loop over points to test

In [ ]:
for (Int_t i = 0; i < parameterScan->numEntries(); ++i) {
//    cout << "on parameter point " << i << " out of " << parameterScan->numEntries() << endl;
// get a parameter point from the list of points to test.
tmpPoint = (RooArgSet *)parameterScan->get(i)->clone("temp");

TMarker *mark = new TMarker(tmpPoint->getRealValue("mu"), 1, 25);
if (interval->IsInInterval(*tmpPoint))
mark->SetMarkerColor(kBlue);
else
mark->SetMarkerColor(kRed);

mark->Draw("s");
// delete tmpPoint;
//    delete mark;
}
t.Stop();
t.Print();


Draw all canvases

In [ ]:
gROOT->GetListOfCanvases()->Draw()