Difference between TFeldmanCousins and StandardHypoTestInvDe

Hi everyone,

I compute the upper limit at 95%CL with

  • the observed number = 1
  • and the expected background number = 0.145

Using TFeldmanCousins, the upper limit is 5, while using StandardHypoTestInvDemo.C in roostats, the upper limit is 2.94 . The value given by TFeldmanCousins is more compatible with the value given in F-C paper (prd.aps.org/abstract/PRD/v57/i7/p3873_1).

To compute this limit by StandardHypoTestInvDemo.C, I do :

.L StandardHypoTestInvDemo.C

StandardHypoTestInvDemo(“PoissonModelWithBackg.root”,“w”,“ModelConfig”,"",“data”,0,2, false, 100, 0, 10, 1000,false, 0)

And the model is constructed as below.

I do not know why there is this difference. What is wrong here?
Thank you very much for your explanations.

Cheers,
Viet Nga

-----------Model of signal and background---------

#include “RooPlot.h”
#include “RooAbsPdf.h”
#include “RooWorkspace.h”
#include “RooDataSet.h”
#include “RooGlobalFunc.h”
#include “RooFitResult.h”
#include “RooRandom.h”
#include “RooAbsReal.h”

#include “RooStats/RooStatsUtils.h”
#include “RooStats/ProfileLikelihoodCalculator.h”
#include “RooStats/LikelihoodInterval.h”
#include “RooStats/LikelihoodIntervalPlot.h”
#include “RooStats/BayesianCalculator.h”
#include “RooStats/MCMCCalculator.h”
#include “RooStats/MCMCInterval.h”
#include “RooStats/MCMCIntervalPlot.h”
#include “RooStats/ProposalHelper.h”
#include “RooStats/SimpleInterval.h”
#include “RooStats/FeldmanCousins.h”
#include “RooStats/PointSetInterval.h”
#include “RooStats/ToyMCSampler.h”
#include "RooStats/ProfileLikelihoodTestStat.h"
using namespace std;
using namespace RooFit;
using namespace RooFit ;

void stau124_noSyst(){
gROOT->Reset();
int nobs = 1; // number of observed events
double b = 0.145; // number of background events

RooWorkspace w(“w”,kTRUE);

// make Poisson model * Gaussian constraint
w.factory(“sum:nexp(s[1,0,50],b[0.145,0,10])”);

// Poisson of (n | s+b)
w.factory(“Poisson:pdf(nobs[0,60],nexp)”);

//===========CHANGE===========================
w.factory(“PROD:model(pdf)”);

//set value of observed events
RooRealVar * obs = w.var(“nobs”);
obs->setVal(nobs);
// make data set with the namber of observed events
RooDataSet data(“data”,"", *obs );
data.add(*obs );
w.import(data);

w.Print();

RooStats::ModelConfig mc(“ModelConfig”,&w);
mc.SetPdf(*w.pdf(“model”));
mc.SetParametersOfInterest(*w.var(“s”));
mc.SetObservables(*w.var(“nobs”));
mc.SetNuisanceParameters(*w.var(“b”));
mc.SetSnapshot(*w.var(“s”));

mc.Print();
// import model in the workspace
w.import(mc);
TString fileName = “PoissonModelWithBackg.root”;

// write workspace in the file (recreate file if already existing)
w.writeToFile(fileName, true);
}

Hi,

Your model description is not correct. You have defined b as b[0.145,0,10]. This makes b not constant but when compute the test statistics it can vary in that range.
Just define the number of expected events as following:

 w.factory("sum:nexp(s[1,0,50],b[0.145])"); 

or call later

w.var("b")->setConstant(true);

After doing this you will get a value much closer to 5. (I got 4.9)

Best Regards

Lorenzo

:smiley:
Thank you Lorenzo.
Cheers,
Viet Nga