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 Sunday, February 16, 2020 at 04:05 AM.

In [1]:
%%cpp -d
#include "RooRealVar.h"
#include "RooProdPdf.h"
#include "RooWorkspace.h"
#include "RooDataSet.h"
#include "TCanvas.h"
#include "TH1.h"
In [2]:
%%cpp -d
// This is a workaround to make sure the namespace is used inside functions
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 [3]:
RooWorkspace *w1 = new RooWorkspace("w", true);
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)");
RooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby 
                Copyright (C) 2000-2013 NIKHEF, University of California & Stanford University
                All rights reserved, please read http://roofit.sourceforge.net/license.txt

[#1] INFO:ObjectHandling -- RooWorkspace::exportToCint(w) INFO: references to all objects in this workspace will be created in CINT in 'namespace w'

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

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

Plot it, blue is averaged model, red is b known exactly

In [5]:
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]_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_NORM[x]]_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 [6]:
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]_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_X_px]_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]_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]_denominator_Int[b]) using numeric integrator RooIntegrator1D to calculate Int(b)
[#1] INFO:NumericIntegration -- RooRealIntegral::init([py_X_prior_b]_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]_X_px_cdf_Int[x_prime|CDF]_Norm[x_prime]]_Int[b|CDF]) using numeric integrator RooIntegrator1D to calculate Int(b)
Hybrid p-value = 0.999226
Z_Gamma Significance  = 3.1655

Analytic z_bi

In [7]:
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 [8]:
%jsroot on
gROOT->GetListOfCanvases()->Draw()