ROOT logo

From $ROOTSYS/tutorials/roostats/rs101_limitexample.C

/////////////////////////////////////////////////////////////////////////
//
// 'Limit Example' RooStats tutorial macro #101
// author: Kyle Cranmer
// date June. 2009
//
// This tutorial shows an example of creating a simple
// model for a number counting experiment with uncertainty
// on both the background rate and signal efficeincy. We then 
// use a Confidence Interval Calculator to set a limit on the signal.
//
//
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif

#include "RooProfileLL.h"
#include "RooAbsPdf.h"
#include "RooStats/HypoTestResult.h"
#include "RooRealVar.h"
#include "RooPlot.h"
#include "RooDataSet.h"
#include "RooTreeDataStore.h"
#include "TTree.h"
#include "TCanvas.h"
#include "TLine.h"
#include "TStopwatch.h"

#include "RooStats/ProfileLikelihoodCalculator.h"
#include "RooStats/MCMCCalculator.h"
#include "RooStats/UniformProposal.h"
#include "RooStats/FeldmanCousins.h"
#include "RooStats/NumberCountingPdfFactory.h"
#include "RooStats/ConfInterval.h"
#include "RooStats/PointSetInterval.h"
#include "RooStats/LikelihoodInterval.h"
#include "RooStats/LikelihoodIntervalPlot.h"
#include "RooStats/RooStatsUtils.h"
#include "RooStats/ModelConfig.h"
#include "RooStats/MCMCInterval.h"
#include "RooStats/MCMCIntervalPlot.h"
#include "RooStats/ProposalFunction.h"
#include "RooStats/ProposalHelper.h"
#include "RooFitResult.h"


// use this order for safety on library loading
using namespace RooFit ;
using namespace RooStats ;


void rs101_limitexample()
{
  /////////////////////////////////////////
  // An example of setting a limit in a number counting experiment with uncertainty on background and signal
  /////////////////////////////////////////

  // to time the macro
  TStopwatch t;
  t.Start();

  /////////////////////////////////////////
  // The Model building stage
  /////////////////////////////////////////
  RooWorkspace* wspace = new RooWorkspace();
  wspace->factory("Poisson::countingModel(obs[150,0,300], sum(s[50,0,120]*ratioSigEff[1.,0,2.],b[100,0,300]*ratioBkgEff[1.,0.,2.]))"); // counting model
  //  wspace->factory("Gaussian::sigConstraint(ratioSigEff,1,0.05)"); // 5% signal efficiency uncertainty
  //  wspace->factory("Gaussian::bkgConstraint(ratioBkgEff,1,0.1)"); // 10% background efficiency uncertainty
  wspace->factory("Gaussian::sigConstraint(1,ratioSigEff,0.05)"); // 5% signal efficiency uncertainty
  wspace->factory("Gaussian::bkgConstraint(1,ratioBkgEff,0.1)"); // 10% background efficiency uncertainty
  wspace->factory("PROD::modelWithConstraints(countingModel,sigConstraint,bkgConstraint)"); // product of terms
  wspace->Print();

  RooAbsPdf* modelWithConstraints = wspace->pdf("modelWithConstraints"); // get the model
  RooRealVar* obs = wspace->var("obs"); // get the observable
  RooRealVar* s = wspace->var("s"); // get the signal we care about
  RooRealVar* b = wspace->var("b"); // get the background and set it to a constant.  Uncertainty included in ratioBkgEff
  b->setConstant();
  RooRealVar* ratioSigEff = wspace->var("ratioSigEff"); // get uncertaint parameter to constrain
  RooRealVar* ratioBkgEff = wspace->var("ratioBkgEff"); // get uncertaint parameter to constrain
  RooArgSet constrainedParams(*ratioSigEff, *ratioBkgEff); // need to constrain these in the fit (should change default behavior)

  // Create an example dataset with 160 observed events
  obs->setVal(160.);
  RooDataSet* data = new RooDataSet("exampleData", "exampleData", RooArgSet(*obs));
  data->add(*obs);

  RooArgSet all(*s, *ratioBkgEff, *ratioSigEff);

  // not necessary
  modelWithConstraints->fitTo(*data, RooFit::Constrain(RooArgSet(*ratioSigEff, *ratioBkgEff)));

  // Now let's make some confidence intervals for s, our parameter of interest
  RooArgSet paramOfInterest(*s);

  ModelConfig modelConfig(new RooWorkspace());
  modelConfig.SetPdf(*modelWithConstraints);
  modelConfig.SetParametersOfInterest(paramOfInterest);


  // First, let's use a Calculator based on the Profile Likelihood Ratio
  //ProfileLikelihoodCalculator plc(*data, *modelWithConstraints, paramOfInterest); 
  ProfileLikelihoodCalculator plc(*data, modelConfig); 
  plc.SetTestSize(.05);
  ConfInterval* lrint = plc.GetInterval();  // that was easy.

  // Let's make a plot
  TCanvas* dataCanvas = new TCanvas("dataCanvas");
  dataCanvas->Divide(2,1);

  dataCanvas->cd(1);
  LikelihoodIntervalPlot plotInt((LikelihoodInterval*)lrint);
  plotInt.SetTitle("Profile Likelihood Ratio and Posterior for S");
  plotInt.Draw();

  // Second, use a Calculator based on the Feldman Cousins technique
  FeldmanCousins fc(*data, modelConfig);
  fc.UseAdaptiveSampling(true);
  fc.FluctuateNumDataEntries(false); // number counting analysis: dataset always has 1 entry with N events observed
  fc.SetNBins(100); // number of points to test per parameter
  fc.SetTestSize(.05);
  //  fc.SaveBeltToFile(true); // optional
  ConfInterval* fcint = NULL;
  fcint = fc.GetInterval();  // that was easy.

  RooFitResult* fit = modelWithConstraints->fitTo(*data, Save(true));

  // Third, use a Calculator based on Markov Chain monte carlo
  // Before configuring the calculator, let's make a ProposalFunction
  // that will achieve a high acceptance rate
  ProposalHelper ph;
  ph.SetVariables((RooArgSet&)fit->floatParsFinal());
  ph.SetCovMatrix(fit->covarianceMatrix());
  ph.SetUpdateProposalParameters(true);
  ph.SetCacheSize(100);
  ProposalFunction* pdfProp = ph.GetProposalFunction();  // that was easy

  MCMCCalculator mc(*data, modelConfig);
  mc.SetNumIters(20000); // steps to propose in the chain
  mc.SetTestSize(.05); // 95% CL
  mc.SetNumBurnInSteps(40); // ignore first N steps in chain as "burn in"
  mc.SetProposalFunction(*pdfProp);
  mc.SetLeftSideTailFraction(0.5);  // find a "central" interval
  MCMCInterval* mcInt = (MCMCInterval*)mc.GetInterval();  // that was easy


  // Get Lower and Upper limits from Profile Calculator
  cout << "Profile lower limit on s = " << ((LikelihoodInterval*) lrint)->LowerLimit(*s) << endl;
  cout << "Profile upper limit on s = " << ((LikelihoodInterval*) lrint)->UpperLimit(*s) << endl;

  // Get Lower and Upper limits from FeldmanCousins with profile construction
  if (fcint != NULL) {
     double fcul = ((PointSetInterval*) fcint)->UpperLimit(*s);
     double fcll = ((PointSetInterval*) fcint)->LowerLimit(*s);
     cout << "FC lower limit on s = " << fcll << endl;
     cout << "FC upper limit on s = " << fcul << endl;
     TLine* fcllLine = new TLine(fcll, 0, fcll, 1);
     TLine* fculLine = new TLine(fcul, 0, fcul, 1);
     fcllLine->SetLineColor(kRed);
     fculLine->SetLineColor(kRed);
     fcllLine->Draw("same");
     fculLine->Draw("same");
     dataCanvas->Update();
  }

  // Plot MCMC interval and print some statistics
  MCMCIntervalPlot mcPlot(*mcInt);
  mcPlot.SetLineColor(kMagenta);
  mcPlot.SetLineWidth(2);
  mcPlot.Draw("same");

  double mcul = mcInt->UpperLimit(*s);
  double mcll = mcInt->LowerLimit(*s);
  cout << "MCMC lower limit on s = " << mcll << endl;
  cout << "MCMC upper limit on s = " << mcul << endl;
  cout << "MCMC Actual confidence level: "
     << mcInt->GetActualConfidenceLevel() << endl;

  // 3-d plot of the parameter points
  dataCanvas->cd(2);
  // also plot the points in the markov chain
  TTree& chain =  ((RooTreeDataStore*) mcInt->GetChainAsDataSet()->store())->tree();
  chain.SetMarkerStyle(6);
  chain.SetMarkerColor(kRed);
  chain.Draw("s:ratioSigEff:ratioBkgEff","weight_MarkovChain_local_","box"); // 3-d box proporional to posterior

  // the points used in the profile construction
  TTree& parameterScan =  ((RooTreeDataStore*) fc.GetPointsToScan()->store())->tree();
  parameterScan.SetMarkerStyle(24);
  parameterScan.Draw("s:ratioSigEff:ratioBkgEff","","same");

  delete wspace;
  delete lrint;
  delete mcInt;
  delete fcint;
  delete data;

  /// print timing info
  t.Stop();
  t.Print();
}
 rs101_limitexample.C:1
 rs101_limitexample.C:2
 rs101_limitexample.C:3
 rs101_limitexample.C:4
 rs101_limitexample.C:5
 rs101_limitexample.C:6
 rs101_limitexample.C:7
 rs101_limitexample.C:8
 rs101_limitexample.C:9
 rs101_limitexample.C:10
 rs101_limitexample.C:11
 rs101_limitexample.C:12
 rs101_limitexample.C:13
 rs101_limitexample.C:14
 rs101_limitexample.C:15
 rs101_limitexample.C:16
 rs101_limitexample.C:17
 rs101_limitexample.C:18
 rs101_limitexample.C:19
 rs101_limitexample.C:20
 rs101_limitexample.C:21
 rs101_limitexample.C:22
 rs101_limitexample.C:23
 rs101_limitexample.C:24
 rs101_limitexample.C:25
 rs101_limitexample.C:26
 rs101_limitexample.C:27
 rs101_limitexample.C:28
 rs101_limitexample.C:29
 rs101_limitexample.C:30
 rs101_limitexample.C:31
 rs101_limitexample.C:32
 rs101_limitexample.C:33
 rs101_limitexample.C:34
 rs101_limitexample.C:35
 rs101_limitexample.C:36
 rs101_limitexample.C:37
 rs101_limitexample.C:38
 rs101_limitexample.C:39
 rs101_limitexample.C:40
 rs101_limitexample.C:41
 rs101_limitexample.C:42
 rs101_limitexample.C:43
 rs101_limitexample.C:44
 rs101_limitexample.C:45
 rs101_limitexample.C:46
 rs101_limitexample.C:47
 rs101_limitexample.C:48
 rs101_limitexample.C:49
 rs101_limitexample.C:50
 rs101_limitexample.C:51
 rs101_limitexample.C:52
 rs101_limitexample.C:53
 rs101_limitexample.C:54
 rs101_limitexample.C:55
 rs101_limitexample.C:56
 rs101_limitexample.C:57
 rs101_limitexample.C:58
 rs101_limitexample.C:59
 rs101_limitexample.C:60
 rs101_limitexample.C:61
 rs101_limitexample.C:62
 rs101_limitexample.C:63
 rs101_limitexample.C:64
 rs101_limitexample.C:65
 rs101_limitexample.C:66
 rs101_limitexample.C:67
 rs101_limitexample.C:68
 rs101_limitexample.C:69
 rs101_limitexample.C:70
 rs101_limitexample.C:71
 rs101_limitexample.C:72
 rs101_limitexample.C:73
 rs101_limitexample.C:74
 rs101_limitexample.C:75
 rs101_limitexample.C:76
 rs101_limitexample.C:77
 rs101_limitexample.C:78
 rs101_limitexample.C:79
 rs101_limitexample.C:80
 rs101_limitexample.C:81
 rs101_limitexample.C:82
 rs101_limitexample.C:83
 rs101_limitexample.C:84
 rs101_limitexample.C:85
 rs101_limitexample.C:86
 rs101_limitexample.C:87
 rs101_limitexample.C:88
 rs101_limitexample.C:89
 rs101_limitexample.C:90
 rs101_limitexample.C:91
 rs101_limitexample.C:92
 rs101_limitexample.C:93
 rs101_limitexample.C:94
 rs101_limitexample.C:95
 rs101_limitexample.C:96
 rs101_limitexample.C:97
 rs101_limitexample.C:98
 rs101_limitexample.C:99
 rs101_limitexample.C:100
 rs101_limitexample.C:101
 rs101_limitexample.C:102
 rs101_limitexample.C:103
 rs101_limitexample.C:104
 rs101_limitexample.C:105
 rs101_limitexample.C:106
 rs101_limitexample.C:107
 rs101_limitexample.C:108
 rs101_limitexample.C:109
 rs101_limitexample.C:110
 rs101_limitexample.C:111
 rs101_limitexample.C:112
 rs101_limitexample.C:113
 rs101_limitexample.C:114
 rs101_limitexample.C:115
 rs101_limitexample.C:116
 rs101_limitexample.C:117
 rs101_limitexample.C:118
 rs101_limitexample.C:119
 rs101_limitexample.C:120
 rs101_limitexample.C:121
 rs101_limitexample.C:122
 rs101_limitexample.C:123
 rs101_limitexample.C:124
 rs101_limitexample.C:125
 rs101_limitexample.C:126
 rs101_limitexample.C:127
 rs101_limitexample.C:128
 rs101_limitexample.C:129
 rs101_limitexample.C:130
 rs101_limitexample.C:131
 rs101_limitexample.C:132
 rs101_limitexample.C:133
 rs101_limitexample.C:134
 rs101_limitexample.C:135
 rs101_limitexample.C:136
 rs101_limitexample.C:137
 rs101_limitexample.C:138
 rs101_limitexample.C:139
 rs101_limitexample.C:140
 rs101_limitexample.C:141
 rs101_limitexample.C:142
 rs101_limitexample.C:143
 rs101_limitexample.C:144
 rs101_limitexample.C:145
 rs101_limitexample.C:146
 rs101_limitexample.C:147
 rs101_limitexample.C:148
 rs101_limitexample.C:149
 rs101_limitexample.C:150
 rs101_limitexample.C:151
 rs101_limitexample.C:152
 rs101_limitexample.C:153
 rs101_limitexample.C:154
 rs101_limitexample.C:155
 rs101_limitexample.C:156
 rs101_limitexample.C:157
 rs101_limitexample.C:158
 rs101_limitexample.C:159
 rs101_limitexample.C:160
 rs101_limitexample.C:161
 rs101_limitexample.C:162
 rs101_limitexample.C:163
 rs101_limitexample.C:164
 rs101_limitexample.C:165
 rs101_limitexample.C:166
 rs101_limitexample.C:167
 rs101_limitexample.C:168
 rs101_limitexample.C:169
 rs101_limitexample.C:170
 rs101_limitexample.C:171
 rs101_limitexample.C:172
 rs101_limitexample.C:173
 rs101_limitexample.C:174
 rs101_limitexample.C:175
 rs101_limitexample.C:176
 rs101_limitexample.C:177
 rs101_limitexample.C:178
 rs101_limitexample.C:179
 rs101_limitexample.C:180
 rs101_limitexample.C:181
 rs101_limitexample.C:182
 rs101_limitexample.C:183
 rs101_limitexample.C:184
 rs101_limitexample.C:185
 rs101_limitexample.C:186
 rs101_limitexample.C:187
 rs101_limitexample.C:188
 rs101_limitexample.C:189
 rs101_limitexample.C:190
 rs101_limitexample.C:191
 rs101_limitexample.C:192
 rs101_limitexample.C:193
 rs101_limitexample.C:194
 rs101_limitexample.C:195
 rs101_limitexample.C:196
 rs101_limitexample.C:197
 rs101_limitexample.C:198
 rs101_limitexample.C:199
 rs101_limitexample.C:200
 rs101_limitexample.C:201
 rs101_limitexample.C:202
 rs101_limitexample.C:203
 rs101_limitexample.C:204