Zbi_Zgamma¶

Demonstrate Z_Bi = Z_Gamma

Author: Kyle Cranmer, Wouter Verkerke
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Monday, May 29, 2023 at 09:33 AM.

In [1]:

%%cpp -d
#include "RooRealVar.h"
#include "RooProdPdf.h"
#include "RooWorkspace.h"
#include "RooDataSet.h"
#include "TCanvas.h"
#include "TH1.h"

using namespace RooFit;
using namespace RooStats;

Make model for prototype on/off problem Pois(x | s+b) * Pois(y | tau b ) for Z_Gamma, use uniform prior on b.

In [2]:

RooWorkspace *w1 = new 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)");

construct the Bayesian-averaged model (eg. a projection pdf) p'(x|s) = \int db p(x|s+b) * [ p(y|b) * prior(b) ]

In [3]:

w1->factory("PROJ::averagedModel(PROD::foo(px|b,py,prior_b),b)");

[#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

In [4]:

RooPlot *frame = w1->var("x")->frame();
w1->pdf("averagedModel")->plotOn(frame);
w1->pdf("px")->plotOn(frame, LineColor(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)

Info in <TCanvas::MakeDefCanvas>:  created default TCanvas with name c1

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

In [5]:

w1->var("y")->setVal(100);
w1->var("x")->setVal(150);
RooAbsReal *cdf = w1->pdf("averagedModel")->createCdf(*w1->var("x"));
cdf->getVal(); // get ugly print messages out of the way

cout << "Hybrid p-value = " << cdf->getVal() << endl;
cout << "Z_Gamma Significance  = " << PValueToSignificance(1 - cdf->getVal()) << endl;

[#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)
Hybrid p-value = 0.999226
Z_Gamma Significance  = 3.1655

analytic Z_Bi

In [6]:

double Z_Bi = NumberCountingUtils::BinomialWithTauObsZ(150, 100, 1);
std::cout << "Z_Bi significance estimation: " << Z_Bi << std::endl;

Z_Bi significance estimation: 3.10804

OUTPUT Hybrid p-value = 0.999058 Z_Gamma Significance = 3.10804 Z_Bi significance estimation: 3.10804

Draw all canvases

In [7]:

%jsroot on
gROOT->GetListOfCanvases()->Draw()