Demonstrate Z_Bi = Z_Gamma
Author: Artem Busorgin, Kyle Cranmer and Wouter Verkerke (C++ version)
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Tuesday, March 19, 2024 at 07:18 PM.
import ROOT
Make model for prototype on/off problem Pois(x | s+b) * Pois(y | tau b ) for Z_Gamma, use uniform prior on b.
w1 = ROOT.RooWorkspace("w")
w1.factory("Poisson::px(x[150,0,500],sum::splusb(s[0,0,100],b[100,0,300]))")
w1.factory("Poisson::py(y[100,0,500],prod::taub(tau[1.],b))")
w1.factory("Uniform::prior_b(b)")
<cppyy.gbl.RooUniform object at 0xaf02380>
construct the Bayesian-averaged model (eg. a projection pdf) p'(x|s) = \int db p(x|s+b) * [ p(y|b) * prior(b) ]
w1.factory("PROJ::averagedModel(PROD::foo(px|b,py,prior_b),b)")
c = ROOT.TCanvas()
[#1] INFO:NumericIntegration -- RooRealIntegral::init([py_X_prior_b_X_px]_Norm[b]_denominator_Int[b]) using numeric integrator RooIntegrator1D to calculate Int(b)
plot it, blue is averaged model, red is b known exactly
frame = w1["x"].frame()
w1["averagedModel"].plotOn(frame)
w1["px"].plotOn(frame, LineColor=ROOT.kRed)
frame.Draw()
[#1] INFO:NumericIntegration -- RooRealIntegral::init([py_X_prior_b]_Norm[b]_denominator_Int[b]) using numeric integrator RooIntegrator1D to calculate Int(b) [#1] INFO:NumericIntegration -- RooRealIntegral::init([[py_X_prior_b]_Norm[b]_X_px_NORM[x]]_Int[b]) using numeric integrator RooIntegrator1D to calculate Int(b) [#1] INFO:NumericIntegration -- RooRealIntegral::init([py_X_prior_b]_Norm[b]_denominator_Int[b]) using numeric integrator RooIntegrator1D to calculate Int(b)
compare analytic calculation of Z_Bi with the numerical RooFit implementation of Z_Gamma for an example with x = 150, y = 100
numeric RooFit Z_Gamma
w1["y"].setVal(100)
w1["x"].setVal(150)
cdf = w1["averagedModel"].createCdf(w1["x"])
cdf.getVal() # get ugly print messages out of the way
print("Hybrid p-value = ", cdf.getVal())
print("Z_Gamma Significance = ", ROOT.RooStats.PValueToSignificance(1 - cdf.getVal()))
Hybrid p-value = 0.9992259057034769 Z_Gamma Significance = 3.165495870670026 [#1] INFO:NumericIntegration -- RooRealIntegral::init([py_X_prior_b]_Norm[b]_denominator_Int[b]) using numeric integrator RooIntegrator1D to calculate Int(b) [#1] INFO:NumericIntegration -- RooRealIntegral::init([py_X_prior_b]_Norm[b]_denominator_Int[b]) using numeric integrator RooIntegrator1D to calculate Int(b) [#1] INFO:NumericIntegration -- RooRealIntegral::init([[py_X_prior_b]_Norm[b]_X_px_cdf_NORM[x_prime]]_Int[b]) using numeric integrator RooIntegrator1D to calculate Int(b) [#1] INFO:NumericIntegration -- RooRealIntegral::init([py_X_prior_b]_Norm[b]_denominator_Int[b]) using numeric integrator RooIntegrator1D to calculate Int(b) [#1] INFO:NumericIntegration -- RooRealIntegral::init([[py_X_prior_b]_Norm[b]_X_px_cdf_Int[x_prime|CDF]_Norm[x_prime]]_Int[b|CDF]) using numeric integrator RooIntegrator1D to calculate Int(b) [#1] INFO:NumericIntegration -- RooRealIntegral::init([py_X_prior_b]_Norm[b]_denominator_Int[b]) using numeric integrator RooIntegrator1D to calculate Int(b)
analytic Z_Bi
Z_Bi = ROOT.RooStats.NumberCountingUtils.BinomialWithTauObsZ(150, 100, 1)
print("Z_Bi significance estimation: ", Z_Bi)
c.SaveAs("Zbi_Zgamma.png")
Z_Bi significance estimation: 3.1080438957471137
Info in <TCanvas::Print>: png file Zbi_Zgamma.png has been created
Draw all canvases
%jsroot on
from ROOT import gROOT
gROOT.GetListOfCanvases().Draw()