ROOT logo

From $ROOTSYS/tutorials/roostats/rs801_HypoTestInverterOriginal.C

/////////////////////////////////////////////////////////////////////////
//
// 'Hypothesis Test Inversion' RooStats tutorial macro #801
// author: Gregory Schott
// date Sep 2009
//
// This tutorial shows an example of using the HypoTestInverterOriginal class 
//
/////////////////////////////////////////////////////////////////////////

#include "RooRealVar.h"
#include "RooConstVar.h"
#include "RooProdPdf.h"
#include "RooWorkspace.h"
#include "RooDataSet.h"
#include "RooPolynomial.h"
#include "RooAddPdf.h"
#include "RooExtendPdf.h"

#include "RooStats/HypoTestInverterOriginal.h"
#include "RooStats/HypoTestInverterResult.h"
#include "RooStats/HypoTestInverterPlot.h"
#include "RooStats/HybridCalculatorOriginal.h"

#include "TGraphErrors.h"

using namespace RooFit;
using namespace RooStats;


void rs801_HypoTestInverterOriginal()
{
  // prepare the model
  RooRealVar lumi("lumi","luminosity",1);
  RooRealVar r("r","cross-section ratio",3.74,0,50);
  RooFormulaVar ns("ns","1*r*lumi",RooArgList(lumi,r));
  RooRealVar nb("nb","background yield",1);
  RooRealVar x("x","dummy observable",0,1);
  RooConstVar p0(RooFit::RooConst(0));
  RooPolynomial flatPdf("flatPdf","flat PDF",x,p0);
  RooAddPdf totPdf("totPdf","S+B model",RooArgList(flatPdf,flatPdf),RooArgList(ns,nb));
  RooExtendPdf bkgPdf("bkgPdf","B-only model",flatPdf,nb);
  RooDataSet* data = totPdf.generate(x,1);

  // prepare the calculator
  HybridCalculatorOriginal myhc(*data, totPdf, bkgPdf,0,0);
  myhc.SetTestStatistic(2);
  myhc.SetNumberOfToys(1000);
  myhc.UseNuisance(false);                            

  // run the hypothesis-test invertion
  HypoTestInverterOriginal myInverter(myhc,r);
  myInverter.SetTestSize(0.10);
  myInverter.UseCLs(true);
  // myInverter.RunFixedScan(5,1,6);
  // scan for a 95% UL
  myInverter.RunAutoScan(3.,5,myInverter.Size()/2,0.005);  
  // run an alternative autoscan algorithm 
  // myInverter.RunAutoScan(1,6,myInverter.Size()/2,0.005,1);  
  //myInverter.RunOnePoint(3.9);


  HypoTestInverterResult* results = myInverter.GetInterval();

  HypoTestInverterPlot myInverterPlot("myInverterPlot","",results);
  TGraphErrors* gr1 = myInverterPlot.MakePlot();
  gr1->Draw("ALP");

  double ulError = results->UpperLimitEstimatedError();

  double upperLimit = results->UpperLimit();
  std::cout << "The computed upper limit is: " << upperLimit << std::endl;
  std::cout << "an estimated error on this upper limit is: " << ulError << std::endl;
  // expected result: 4.10
}
int main() { 
   rs801_HypoTestInverter();
}
 rs801_HypoTestInverterOriginal.C:1
 rs801_HypoTestInverterOriginal.C:2
 rs801_HypoTestInverterOriginal.C:3
 rs801_HypoTestInverterOriginal.C:4
 rs801_HypoTestInverterOriginal.C:5
 rs801_HypoTestInverterOriginal.C:6
 rs801_HypoTestInverterOriginal.C:7
 rs801_HypoTestInverterOriginal.C:8
 rs801_HypoTestInverterOriginal.C:9
 rs801_HypoTestInverterOriginal.C:10
 rs801_HypoTestInverterOriginal.C:11
 rs801_HypoTestInverterOriginal.C:12
 rs801_HypoTestInverterOriginal.C:13
 rs801_HypoTestInverterOriginal.C:14
 rs801_HypoTestInverterOriginal.C:15
 rs801_HypoTestInverterOriginal.C:16
 rs801_HypoTestInverterOriginal.C:17
 rs801_HypoTestInverterOriginal.C:18
 rs801_HypoTestInverterOriginal.C:19
 rs801_HypoTestInverterOriginal.C:20
 rs801_HypoTestInverterOriginal.C:21
 rs801_HypoTestInverterOriginal.C:22
 rs801_HypoTestInverterOriginal.C:23
 rs801_HypoTestInverterOriginal.C:24
 rs801_HypoTestInverterOriginal.C:25
 rs801_HypoTestInverterOriginal.C:26
 rs801_HypoTestInverterOriginal.C:27
 rs801_HypoTestInverterOriginal.C:28
 rs801_HypoTestInverterOriginal.C:29
 rs801_HypoTestInverterOriginal.C:30
 rs801_HypoTestInverterOriginal.C:31
 rs801_HypoTestInverterOriginal.C:32
 rs801_HypoTestInverterOriginal.C:33
 rs801_HypoTestInverterOriginal.C:34
 rs801_HypoTestInverterOriginal.C:35
 rs801_HypoTestInverterOriginal.C:36
 rs801_HypoTestInverterOriginal.C:37
 rs801_HypoTestInverterOriginal.C:38
 rs801_HypoTestInverterOriginal.C:39
 rs801_HypoTestInverterOriginal.C:40
 rs801_HypoTestInverterOriginal.C:41
 rs801_HypoTestInverterOriginal.C:42
 rs801_HypoTestInverterOriginal.C:43
 rs801_HypoTestInverterOriginal.C:44
 rs801_HypoTestInverterOriginal.C:45
 rs801_HypoTestInverterOriginal.C:46
 rs801_HypoTestInverterOriginal.C:47
 rs801_HypoTestInverterOriginal.C:48
 rs801_HypoTestInverterOriginal.C:49
 rs801_HypoTestInverterOriginal.C:50
 rs801_HypoTestInverterOriginal.C:51
 rs801_HypoTestInverterOriginal.C:52
 rs801_HypoTestInverterOriginal.C:53
 rs801_HypoTestInverterOriginal.C:54
 rs801_HypoTestInverterOriginal.C:55
 rs801_HypoTestInverterOriginal.C:56
 rs801_HypoTestInverterOriginal.C:57
 rs801_HypoTestInverterOriginal.C:58
 rs801_HypoTestInverterOriginal.C:59
 rs801_HypoTestInverterOriginal.C:60
 rs801_HypoTestInverterOriginal.C:61
 rs801_HypoTestInverterOriginal.C:62
 rs801_HypoTestInverterOriginal.C:63
 rs801_HypoTestInverterOriginal.C:64
 rs801_HypoTestInverterOriginal.C:65
 rs801_HypoTestInverterOriginal.C:66
 rs801_HypoTestInverterOriginal.C:67
 rs801_HypoTestInverterOriginal.C:68
 rs801_HypoTestInverterOriginal.C:69
 rs801_HypoTestInverterOriginal.C:70
 rs801_HypoTestInverterOriginal.C:71
 rs801_HypoTestInverterOriginal.C:72
 rs801_HypoTestInverterOriginal.C:73
 rs801_HypoTestInverterOriginal.C:74
 rs801_HypoTestInverterOriginal.C:75
 rs801_HypoTestInverterOriginal.C:76
 rs801_HypoTestInverterOriginal.C:77
 rs801_HypoTestInverterOriginal.C:78
 rs801_HypoTestInverterOriginal.C:79