Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
rs401c_FeldmanCousins.py
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## \date July 2022
17## \authors Artem Busorgin, Kyle Cranmer (C++ version)
18
19import ROOT
20
21# to time the macro... about 30 s
22t = ROOT.TStopwatch()
23t.Start()
24
25# make a simple model
26x = ROOT.RooRealVar("x", "", 1, 0, 50)
27mu = ROOT.RooRealVar("mu", "", 2.5, 0, 15) # with a limit on mu>=0
28b = ROOT.RooConstVar("b", "", 3.0)
29mean = ROOT.RooAddition("mean", "", [mu, b])
30pois = ROOT.RooPoisson("pois", "", x, mean)
31parameters = {mu}
32
33# create a toy dataset
34data = pois.generate({x}, 1)
35data.Print("v")
36
37dataCanvas = ROOT.TCanvas("dataCanvas")
38frame = x.frame()
39data.plotOn(frame)
40frame.Draw()
41dataCanvas.Update()
42
43w = ROOT.RooWorkspace()
44modelConfig = ROOT.RooStats.ModelConfig("poissonProblem", w)
45modelConfig.SetPdf(pois)
46modelConfig.SetParametersOfInterest(parameters)
47modelConfig.SetObservables({x})
48w.Print()
49
50# show use of Feldman-Cousins
51fc = ROOT.RooStats.FeldmanCousins(data, modelConfig)
52fc.SetTestSize(0.05) # set size of test
53fc.UseAdaptiveSampling(True)
54fc.FluctuateNumDataEntries(False) # number counting analysis: dataset always has 1 entry with N events observed
55fc.SetNBins(100) # number of points to test per parameter
56
57# use the Feldman-Cousins tool
58interval = fc.GetInterval()
59
60# make a canvas for plots
61intervalCanvas = ROOT.TCanvas("intervalCanvas")
62
63print("is this point in the interval? ", interval.IsInInterval(parameters))
64print("interval is [{}, {}]".format(interval.LowerLimit(mu), interval.UpperLimit(mu)))
65
66# using 200 bins it takes 1 min and the answer is
67# interval is [0.2625, 10.6125] with a step size of .075
68# The interval in Feldman & Cousins's original paper is [.29, 10.81]
69# Phys.Rev.D57:3873-3889,1998.
70
71# No dedicated plotting class yet, so do it by hand:
72parameterScan = fc.GetPointsToScan()
73hist = parameterScan.createHistogram("mu", ROOT.RooFit.Binning(30))
74hist.Draw()
75
76marks = []
77# loop over points to test
78for i in range(parameterScan.numEntries()):
79 # get a parameter point from the list of points to test.
80 tmpPoint = parameterScan.get(i).clone("temp")
81
82 mark = ROOT.TMarker(tmpPoint.getRealValue("mu"), 1, 25)
83 if interval.IsInInterval(tmpPoint):
84 mark.SetMarkerColor(ROOT.kBlue)
85 else:
86 mark.SetMarkerColor(ROOT.kRed)
87
88 mark.Draw("s")
89 marks.append(mark)
90
91t.Stop()
92t.Print()
93
94dataCanvas.SaveAs("rs401c_FeldmanCousins_data.png")
95intervalCanvas.SaveAs("rs401c_FeldmanCousins_hist.png")
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t format