rs401c_FeldmanCousins.C File Reference

Detailed Description

Produces an interval on the mean signal in a number counting experiment with known background using the Feldman-Cousins technique.

Using the RooStats FeldmanCousins tool with 200 bins it takes 1 min and the interval is [0.2625, 10.6125] with a step size of 0.075. The interval in Feldman & Cousins's original paper is [.29, 10.81] Phys.Rev.D57:3873-3889,1998.

Processing /mnt/build/workspace/root-makedoc-v612/rootspi/rdoc/src/v6-12-00-patches/tutorials/roostats/rs401c_FeldmanCousins.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

DataStore poisData (Generated From )
  Contains 1 entries
  Observables: 
    1)  x = 7  L(0 - 50)  ""

RooWorkspace()  contents

variables
---------
(mu,x)

p.d.f.s
-------
RooPoisson::pois[ x=x mean=mean ] = 0.0224772

functions
--------
RooAddition::mean[ mu + b ] = 5.5

named sets
----------
poissonProblem_Observables:(x)
poissonProblem_POI:(mu)


=== Using the following for poissonProblem ===
Observables:             RooArgSet:: = (x)
Parameters of Interest:  RooArgSet:: = (mu)
PDF:                     RooPoisson::pois[ x=x mean=mean ] = 0.0224772

FeldmanCousins: ntoys per point: adaptive
FeldmanCousins: nEvents per toy will not fluctuate, will always be 1
FeldmanCousins: Model has no nuisance parameters
FeldmanCousins: # points to test = 100
NeymanConstruction: Prog: 1/100 total MC = 240 this test stat = 1.83324
 mu=0.075 [-1e+30, 1.08573]  in interval = 0

NeymanConstruction: Prog: 2/100 total MC = 80 this test stat = 1.64984
 mu=0.225 [-1e+30, 0.949959]  in interval = 0

NeymanConstruction: Prog: 3/100 total MC = 80 this test stat = 1.4816
 mu=0.375 [-1e+30, 0.827185]  in interval = 0

NeymanConstruction: Prog: 4/100 total MC = 240 this test stat = 1.32721
 mu=0.525 [-1e+30, 1.32721]  in interval = 1

NeymanConstruction: Prog: 5/100 total MC = 80 this test stat = 1.1855
 mu=0.675 [-1e+30, 2.73604]  in interval = 1

NeymanConstruction: Prog: 6/100 total MC = 240 this test stat = 1.05546
 mu=0.825 [-1e+30, 1.72806]  in interval = 1

NeymanConstruction: Prog: 7/100 total MC = 80 this test stat = 0.936198
 mu=0.975 [-1e+30, 2.32979]  in interval = 1

NeymanConstruction: Prog: 8/100 total MC = 80 this test stat = 0.826909
 mu=1.125 [-1e+30, 1.424]  in interval = 1

NeymanConstruction: Prog: 9/100 total MC = 80 this test stat = 0.726882
 mu=1.275 [-1e+30, 1.28826]  in interval = 1

NeymanConstruction: Prog: 10/100 total MC = 80 this test stat = 0.635479
 mu=1.425 [-1e+30, 1.81459]  in interval = 1

NeymanConstruction: Prog: 11/100 total MC = 80 this test stat = 0.552124
 mu=1.575 [-1e+30, 1.66456]  in interval = 1

NeymanConstruction: Prog: 12/100 total MC = 80 this test stat = 0.476298
 mu=1.725 [-1e+30, 1.725]  in interval = 1

NeymanConstruction: Prog: 13/100 total MC = 80 this test stat = 0.40753
 mu=1.875 [-1e+30, 2.05965]  in interval = 1

NeymanConstruction: Prog: 14/100 total MC = 80 this test stat = 0.345393
 mu=2.025 [-1e+30, 1.9066]  in interval = 1

NeymanConstruction: Prog: 15/100 total MC = 80 this test stat = 0.289496
 mu=2.175 [-1e+30, 1.76244]  in interval = 1

NeymanConstruction: Prog: 16/100 total MC = 80 this test stat = 0.239482
 mu=2.325 [-1e+30, 1.7512]  in interval = 1

NeymanConstruction: Prog: 17/100 total MC = 80 this test stat = 0.195006
 mu=2.475 [-1e+30, 1.27184]  in interval = 1

NeymanConstruction: Prog: 18/100 total MC = 80 this test stat = 0.155824
 mu=2.625 [-1e+30, 1.99639]  in interval = 1

NeymanConstruction: Prog: 19/100 total MC = 80 this test stat = 0.121603
 mu=2.775 [-1e+30, 1.86293]  in interval = 1

NeymanConstruction: Prog: 20/100 total MC = 80 this test stat = 0.0921062
 mu=2.925 [-1e+30, 1.73086]  in interval = 1

NeymanConstruction: Prog: 21/100 total MC = 80 this test stat = 0.0670922
 mu=3.075 [-1e+30, 1.60585]  in interval = 1

NeymanConstruction: Prog: 22/100 total MC = 80 this test stat = 0.0463567
 mu=3.225 [-1e+30, 2.10098]  in interval = 1

NeymanConstruction: Prog: 23/100 total MC = 80 this test stat = 0.0296823
 mu=3.375 [-1e+30, 1.96524]  in interval = 1

NeymanConstruction: Prog: 24/100 total MC = 80 this test stat = 0.0168841
 mu=3.525 [-1e+30, 1.97094]  in interval = 1

NeymanConstruction: Prog: 25/100 total MC = 80 this test stat = 0.00778646
 mu=3.675 [-1e+30, 2.07549]  in interval = 1

NeymanConstruction: Prog: 26/100 total MC = 80 this test stat = 0.00222463
 mu=3.825 [-1e+30, 1.59677]  in interval = 1

NeymanConstruction: Prog: 27/100 total MC = 80 this test stat = 4.47494e-05
 mu=3.975 [-1e+30, 2.06902]  in interval = 1

NeymanConstruction: Prog: 28/100 total MC = 80 this test stat = 0.00110295
 mu=4.125 [-1e+30, 2.39501]  in interval = 1

NeymanConstruction: Prog: 29/100 total MC = 80 this test stat = 0.00526396
 mu=4.275 [-1e+30, 1.6175]  in interval = 1

NeymanConstruction: Prog: 30/100 total MC = 80 this test stat = 0.0123995
 mu=4.425 [-1e+30, 1.70627]  in interval = 1

NeymanConstruction: Prog: 31/100 total MC = 80 this test stat = 0.0223862
 mu=4.575 [-1e+30, 1.79627]  in interval = 1

NeymanConstruction: Prog: 32/100 total MC = 80 this test stat = 0.0351383
 mu=4.725 [-1e+30, 2.04934]  in interval = 1

NeymanConstruction: Prog: 33/100 total MC = 80 this test stat = 0.0505187
 mu=4.875 [-1e+30, 2.94484]  in interval = 1

NeymanConstruction: Prog: 34/100 total MC = 80 this test stat = 0.0684391
 mu=5.025 [-1e+30, 3.0571]  in interval = 1

NeymanConstruction: Prog: 35/100 total MC = 80 this test stat = 0.0888053
 mu=5.175 [-1e+30, 2.27954]  in interval = 1

NeymanConstruction: Prog: 36/100 total MC = 80 this test stat = 0.11153
 mu=5.325 [-1e+30, 2.26305]  in interval = 1

NeymanConstruction: Prog: 37/100 total MC = 80 this test stat = 0.136526
 mu=5.475 [-1e+30, 2.03894]  in interval = 1

NeymanConstruction: Prog: 38/100 total MC = 80 this test stat = 0.163716
 mu=5.625 [-1e+30, 1.55152]  in interval = 1

NeymanConstruction: Prog: 39/100 total MC = 80 this test stat = 0.193024
 mu=5.775 [-1e+30, 1.81715]  in interval = 1

NeymanConstruction: Prog: 40/100 total MC = 80 this test stat = 0.224377
 mu=5.925 [-1e+30, 1.71475]  in interval = 1

NeymanConstruction: Prog: 41/100 total MC = 80 this test stat = 0.257707
 mu=6.075 [-1e+30, 2.75427]  in interval = 1

NeymanConstruction: Prog: 42/100 total MC = 80 this test stat = 0.292951
 mu=6.225 [-1e+30, 2.03573]  in interval = 1

NeymanConstruction: Prog: 43/100 total MC = 80 this test stat = 0.330045
 mu=6.375 [-1e+30, 1.92767]  in interval = 1

NeymanConstruction: Prog: 44/100 total MC = 80 this test stat = 0.368932
 mu=6.525 [-1e+30, 2.0545]  in interval = 1

NeymanConstruction: Prog: 45/100 total MC = 80 this test stat = 0.409554
 mu=6.675 [-1e+30, 2.142]  in interval = 1

NeymanConstruction: Prog: 46/100 total MC = 80 this test stat = 0.45186
 mu=6.825 [-1e+30, 2.72294]  in interval = 1

NeymanConstruction: Prog: 47/100 total MC = 80 this test stat = 0.495797
 mu=6.975 [-1e+30, 2.03823]  in interval = 1

NeymanConstruction: Prog: 48/100 total MC = 80 this test stat = 0.541318
 mu=7.125 [-1e+30, 2.41015]  in interval = 1

NeymanConstruction: Prog: 49/100 total MC = 80 this test stat = 0.588375
 mu=7.275 [-1e+30, 1.83449]  in interval = 1

NeymanConstruction: Prog: 50/100 total MC = 80 this test stat = 0.636924
 mu=7.425 [-1e+30, 2.80213]  in interval = 1

NeymanConstruction: Prog: 51/100 total MC = 80 this test stat = 0.686922
 mu=7.575 [-1e+30, 1.82973]  in interval = 1

NeymanConstruction: Prog: 52/100 total MC = 80 this test stat = 0.738329
 mu=7.725 [-1e+30, 1.9093]  in interval = 1

NeymanConstruction: Prog: 53/100 total MC = 80 this test stat = 0.791105
 mu=7.875 [-1e+30, 1.98986]  in interval = 1

NeymanConstruction: Prog: 54/100 total MC = 80 this test stat = 0.845213
 mu=8.025 [-1e+30, 2.81879]  in interval = 1

NeymanConstruction: Prog: 55/100 total MC = 80 this test stat = 0.900613
 mu=8.175 [-1e+30, 2.23216]  in interval = 1

NeymanConstruction: Prog: 56/100 total MC = 80 this test stat = 0.957282
 mu=8.325 [-1e+30, 1.51348]  in interval = 1

NeymanConstruction: Prog: 57/100 total MC = 80 this test stat = 1.01518
 mu=8.475 [-1e+30, 2.32134]  in interval = 1

NeymanConstruction: Prog: 58/100 total MC = 80 this test stat = 1.07427
 mu=8.625 [-1e+30, 4.56136]  in interval = 1

NeymanConstruction: Prog: 59/100 total MC = 80 this test stat = 1.13452
 mu=8.775 [-1e+30, 1.83847]  in interval = 1

NeymanConstruction: Prog: 60/100 total MC = 80 this test stat = 1.19591
 mu=8.925 [-1e+30, 1.80373]  in interval = 1

NeymanConstruction: Prog: 61/100 total MC = 80 this test stat = 1.25841
 mu=9.075 [-1e+30, 2.45893]  in interval = 1

NeymanConstruction: Prog: 62/100 total MC = 80 this test stat = 1.32199
 mu=9.225 [-1e+30, 1.95466]  in interval = 1

NeymanConstruction: Prog: 63/100 total MC = 80 this test stat = 1.38662
 mu=9.375 [-1e+30, 2.24356]  in interval = 1

NeymanConstruction: Prog: 64/100 total MC = 80 this test stat = 1.45228
 mu=9.525 [-1e+30, 2.50319]  in interval = 1

NeymanConstruction: Prog: 65/100 total MC = 80 this test stat = 1.51895
 mu=9.675 [-1e+30, 3.02403]  in interval = 1

NeymanConstruction: Prog: 66/100 total MC = 240 this test stat = 1.5866
 mu=9.825 [-1e+30, 1.60451]  in interval = 1

NeymanConstruction: Prog: 67/100 total MC = 80 this test stat = 1.6552
 mu=9.975 [-1e+30, 3.20706]  in interval = 1

NeymanConstruction: Prog: 68/100 total MC = 80 this test stat = 1.72474
 mu=10.125 [-1e+30, 2.42844]  in interval = 1

[#0] PROGRESS:Generation -- generated toys: 500 / 1440
[#0] PROGRESS:Generation -- generated toys: 1000 / 1440
NeymanConstruction: Prog: 69/100 total MC = 2160 this test stat = 1.79519
 mu=10.275 [-1e+30, 1.79519]  in interval = 1

NeymanConstruction: Prog: 70/100 total MC = 240 this test stat = 1.86654
 mu=10.425 [-1e+30, 1.86654]  in interval = 1

NeymanConstruction: Prog: 71/100 total MC = 80 this test stat = 1.93876
 mu=10.575 [-1e+30, 2.67618]  in interval = 1

NeymanConstruction: Prog: 72/100 total MC = 80 this test stat = 2.01184
 mu=10.725 [-1e+30, 1.40678]  in interval = 0

NeymanConstruction: Prog: 73/100 total MC = 80 this test stat = 2.08575
 mu=10.875 [-1e+30, 1.86151]  in interval = 0

NeymanConstruction: Prog: 74/100 total MC = 240 this test stat = 2.16048
 mu=11.025 [-1e+30, 1.7642]  in interval = 0

NeymanConstruction: Prog: 75/100 total MC = 80 this test stat = 2.23601
 mu=11.175 [-1e+30, 1.59869]  in interval = 0

NeymanConstruction: Prog: 76/100 total MC = 240 this test stat = 2.31232
 mu=11.325 [-1e+30, 1.66448]  in interval = 0

NeymanConstruction: Prog: 77/100 total MC = 80 this test stat = 2.38941
 mu=11.475 [-1e+30, 1.26987]  in interval = 0

NeymanConstruction: Prog: 78/100 total MC = 240 this test stat = 2.46724
 mu=11.625 [-1e+30, 1.40071]  in interval = 0

NeymanConstruction: Prog: 79/100 total MC = 80 this test stat = 2.54582
 mu=11.775 [-1e+30, 1.86704]  in interval = 0

NeymanConstruction: Prog: 80/100 total MC = 80 this test stat = 2.62511
 mu=11.925 [-1e+30, 1.93623]  in interval = 0

NeymanConstruction: Prog: 81/100 total MC = 80 this test stat = 2.70511
 mu=12.075 [-1e+30, 1.43268]  in interval = 0

NeymanConstruction: Prog: 82/100 total MC = 240 this test stat = 2.7858
 mu=12.225 [-1e+30, 1.49357]  in interval = 0

NeymanConstruction: Prog: 83/100 total MC = 80 this test stat = 2.86717
 mu=12.375 [-1e+30, 1.55534]  in interval = 0

NeymanConstruction: Prog: 84/100 total MC = 80 this test stat = 2.94921
 mu=12.525 [-1e+30, 1.12633]  in interval = 0

NeymanConstruction: Prog: 85/100 total MC = 80 this test stat = 3.0319
 mu=12.675 [-1e+30, 2.29398]  in interval = 0

NeymanConstruction: Prog: 86/100 total MC = 80 this test stat = 3.11523
 mu=12.825 [-1e+30, 1.7457]  in interval = 0

NeymanConstruction: Prog: 87/100 total MC = 80 this test stat = 3.1992
 mu=12.975 [-1e+30, 1.8108]  in interval = 0

NeymanConstruction: Prog: 88/100 total MC = 80 this test stat = 3.28378
 mu=13.125 [-1e+30, 2.51757]  in interval = 0

NeymanConstruction: Prog: 89/100 total MC = 240 this test stat = 3.36896
 mu=13.275 [-1e+30, 1.40455]  in interval = 0

NeymanConstruction: Prog: 90/100 total MC = 80 this test stat = 3.45474
 mu=13.425 [-1e+30, 2.01076]  in interval = 0

NeymanConstruction: Prog: 91/100 total MC = 240 this test stat = 3.5411
 mu=13.575 [-1e+30, 2.07896]  in interval = 0

NeymanConstruction: Prog: 92/100 total MC = 80 this test stat = 3.62804
 mu=13.725 [-1e+30, 2.14788]  in interval = 0

NeymanConstruction: Prog: 93/100 total MC = 80 this test stat = 3.71554
 mu=13.875 [-1e+30, 1.16768]  in interval = 0

NeymanConstruction: Prog: 94/100 total MC = 80 this test stat = 3.80359
 mu=14.025 [-1e+30, 1.70402]  in interval = 0

NeymanConstruction: Prog: 95/100 total MC = 240 this test stat = 3.89219
 mu=14.175 [-1e+30, 1.7663]  in interval = 0

NeymanConstruction: Prog: 96/100 total MC = 240 this test stat = 3.98132
 mu=14.325 [-1e+30, 1.32819]  in interval = 0

NeymanConstruction: Prog: 97/100 total MC = 80 this test stat = 4.07097
 mu=14.475 [-1e+30, 1.38336]  in interval = 0

NeymanConstruction: Prog: 98/100 total MC = 80 this test stat = 4.16114
 mu=14.625 [-1e+30, 1.43935]  in interval = 0

NeymanConstruction: Prog: 99/100 total MC = 240 this test stat = 4.25182
 mu=14.775 [-1e+30, 1.49612]  in interval = 0

NeymanConstruction: Prog: 100/100 total MC = 240 this test stat = 4.343
 mu=14.925 [-1e+30, 2.08888]  in interval = 0

[#1] INFO:Eval -- 68 points in interval
is this point in the interval? 0
interval is [0.525, 10.575]
Real time 0:00:03, CP time 3.920

#include "RooGlobalFunc.h"
#include "RooStats/ConfInterval.h"
#include "RooStats/PointSetInterval.h"
#include "RooStats/ConfidenceBelt.h"
#include "RooStats/FeldmanCousins.h"
#include "RooStats/ModelConfig.h"

#include "RooWorkspace.h"
#include "RooDataSet.h"
#include "RooRealVar.h"
#include "RooConstVar.h"
#include "RooAddition.h"

#include "RooDataHist.h"

#include "RooPoisson.h"
#include "RooPlot.h"

#include "TCanvas.h"
#include "TTree.h"
#include "TH1F.h"
#include "TMarker.h"
#include "TStopwatch.h"

#include <iostream>

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


void rs401c_FeldmanCousins()
{
   // to time the macro... about 30 s
   TStopwatch t;
   t.Start();

   // make a simple model
   RooRealVar x("x","", 1,0,50);
   RooRealVar mu("mu","", 2.5,0, 15); // with a limit on mu>=0
   RooConstVar b("b","", 3.);
   RooAddition mean("mean","",RooArgList(mu,b));
   RooPoisson pois("pois", "", x, mean);
   RooArgSet parameters(mu);

   // create a toy dataset
   RooDataSet* data = pois.generate(RooArgSet(x), 1);
   data->Print("v");

   TCanvas* dataCanvas = new TCanvas("dataCanvas");
   RooPlot* frame = x.frame();
   data->plotOn(frame);
   frame->Draw();
   dataCanvas->Update();

   RooWorkspace* w = new RooWorkspace();
   ModelConfig modelConfig("poissonProblem",w);
   modelConfig.SetPdf(pois);
   modelConfig.SetParametersOfInterest(parameters);
   modelConfig.SetObservables(RooArgSet(x));
   w->Print();

   //////// show use of Feldman-Cousins
   RooStats::FeldmanCousins fc(*data,modelConfig);
   fc.SetTestSize(.05); // set size of test
   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

   // use the Feldman-Cousins tool
   PointSetInterval* interval = (PointSetInterval*)fc.GetInterval();

   // make a canvas for plots
   TCanvas* intervalCanvas =  new TCanvas("intervalCanvas");

   std::cout << "is this point in the interval? " <<
      interval->IsInInterval(parameters) << std::endl;

   std::cout << "interval is ["<<
      interval->LowerLimit(mu)  << ", "  <<
      interval->UpperLimit(mu) << "]" << endl;

   // using 200 bins it takes 1 min and the answer is
   // interval is [0.2625, 10.6125] with a step size of .075
   // The interval in Feldman & Cousins's original paper is [.29, 10.81]
   //  Phys.Rev.D57:3873-3889,1998.

   // No dedicated plotting class yet, so do it by hand:

   RooDataHist* parameterScan = (RooDataHist*) fc.GetPointsToScan();
   TH1F* hist = (TH1F*) parameterScan->createHistogram("mu",30);
   hist->Draw();


   RooArgSet* tmpPoint;
   // loop over points to test
   for(Int_t i=0; i<parameterScan->numEntries(); ++i){
      //    cout << "on parameter point " << i << " out of " << parameterScan->numEntries() << endl;
      // get a parameter point from the list of points to test.
      tmpPoint = (RooArgSet*) parameterScan->get(i)->clone("temp");

      TMarker* mark = new TMarker(tmpPoint->getRealValue("mu"), 1, 25);
      if (interval->IsInInterval( *tmpPoint ) )
         mark->SetMarkerColor(kBlue);
      else
         mark->SetMarkerColor(kRed);

      mark->Draw("s");
      //delete tmpPoint;
      //    delete mark;
   }
   t.Stop();
   t.Print();


}

Author

Kyle Cranmer

Definition in file rs401c_FeldmanCousins.C.