Zbi_Zgamma.C File Reference

Detailed Description

Demonstrate Z_Bi = Z_Gamma.

Processing /builddir/build/BUILD/root-6.10.00/tutorials/roostats/Zbi_Zgamma.C...

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'
[#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)
[#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.999058
Z_Gamma Significance  = 3.10804
Z_Bi significance estimation: 3.10804

#include "RooRealVar.h"
#include "RooProdPdf.h"
#include "RooWorkspace.h"
#include "RooDataSet.h"
#include "TCanvas.h"
#include "TH1.h"

using namespace RooFit;
using namespace RooStats;

void Zbi_Zgamma() {

   // Make model for prototype on/off problem
   // Pois(x | s+b) * Pois(y | tau b )
   // for Z_Gamma, use uniform prior on b.
   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)");

   // 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)") ;

   // plot it, blue is averaged model, red is b known exactly
   RooPlot* frame = w1->var("x")->frame() ;
   w1->pdf("averagedModel")->plotOn(frame) ;
   w1->pdf("px")->plotOn(frame,LineColor(kRed)) ;
   frame->Draw() ;

   // 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->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;

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

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

Author

Kyle Cranmer & Wouter Verkerke

Definition in file Zbi_Zgamma.C.