Hi Debasish,
Thanks very much. I think i need some completely different procedure. Ok.
With regards,
Sabya
Hi Debashis!
I am using the macro StandardHypoTestInvDemo.C to calculate the FC limits But I can’t get the same results of TFeldmanCousin.C
The Model Config is below:
using namespace RooFit;
using namespace RooStats;
void CountingModel( double nobs = 1, // number of observed events
double b = 1, // number of background events
double sigmab = 0.001) // relative uncertainty in b
{
RooWorkspace w("w");
// make Poisson model * Gaussian constraint
w.factory("sum:nexp(s[3,0,15],b[1,0,10])");
// Poisson of (n | s+b)
w.factory("Poisson:pdf(nobs[0,50],nexp)");
w.factory("Gaussian:constraint(b0[0,10],b,sigmab[0.01])");
w.factory("PROD:model(pdf,constraint)");
w.var("b")->setVal(b);
w.var("nobs")->setVal(nobs);
w.var("s")->setVal(std::max(1.0,nobs-b));
w.var("b0")->setVal(b);
w.var("b0")->setConstant(true); // needed for being treated as global observables
w.var("sigmab")->setVal(sigmab*b);
ModelConfig mc("ModelConfig",&w);
mc.SetPdf(*w.pdf("model"));
mc.SetParametersOfInterest(*w.var("s"));
mc.SetObservables(*w.var("nobs"));
mc.SetNuisanceParameters(*w.var("b"));
// these are needed for the hypothesis tests
mc.SetSnapshot(*w.var("s"));
mc.SetGlobalObservables(*w.var("b0"));
mc.Print();
// import model in the workspace
w.import(mc);
// make data set with the namber of observed events
RooDataSet data("data","", *w.var("nobs"));
// w.var("nobs")->setVal(3);
data.add(*w.var("nobs") );
// import data set in workspace and save it in a file
w.import(data);
w.Print();
TString fileName = "CountingModel.root";
// write workspace in the file (recreate file if already existing)
w.writeToFile(fileName, true);
}
and the run the below command:
StandardHypoTestInvDemo("CountingModel.root", "w", "ModelConfig", "", "data", 0, 2, false, 20, 0, 6, 1000,true )
But the 95% UL is 3.21. below is the output
But the TFeldmanCousins gives the results [0, 4.14]
So Could you help me? Is there something wrong with my code?
Thank you very much
Hi,
I believe TFeldmanCousin.C code doesn’t include the uncertainty? But StandardHypoTestInvDemo.C takes the uncertainty. So one can use the POLE program and StandardHypoTestInvDemo.C for a better comparison. And to verify the code, maybe someone from ROOT can help if you are giving correctly the arguments to StandardHypoTestInvDemo.C.
Best
Debashis
Hi !
Yes, TFeldmanCousin.C doesn’t include the uncertainty. But if you make the uncertainty smaller and smaller the result should degenerate to the TFeldmanCousin.C.
In fact,No matter how much I change the magnitude of the uncertainty. the results
of StandardHypoTestInvDemo.C is different from the TFeldmanCousin.C.
Hi!
I see that you have used the StandardHypoTestInvDemo.C and get the reasonable results. So could you give me some example of your ModelConfig code ? I want to try and check my codes.
Hi,
It has been some time, I used these codes. You can go over this post also. There were some useful discussions. Also, some working codes are attached in this post.
Hi!
Thanks, I have gone over the posts rencently but I can’t repeat your results in the posts. Is there some updates of StandardHypoTestInvDemo.C in recent years?
Hi,
Can you attach which code you have used? What was my result and what now you are getting? And also how you are using this codes.
I don’t have any idea if StandardHypoTestInvDemo.C is updated in recent years.
Hi !
n_obs = 0, n_bkg = 1.0 sigmab = 0.001
I get the result below:
Time to perform limit scan
Real time 0:00:10, CP time 10.180
The computed upper limit is: 1.51515 +/- 0.12466
Expected upper limits, using the B (alternate) model :
expected limit (median) 4.05556
expected limit (-1 sig) 1.6299
expected limit (+1 sig) 5.84593
expected limit (-2 sig) 0
expected limit (+2 sig) 7.54205
Info in <StandardHypoTestInvDemo>: HypoTestInverterResult has been written in the file Freq_Cls+b_grid_ts2_CountingModel.root
The code is
using namespace RooFit;
using namespace RooStats;
void CountingModel( int nobs = 0, // number of observed events
double b = 1, // number of background events
double sigmab = 0.001) // relative uncertainty in b
{
RooWorkspace w("w");
// make Poisson model * Gaussian constraint
w.factory("sum:nexp(s[3,0,15],b[1,0,10])");
// Poisson of (n | s+b)
w.factory("Poisson:pdf(nobs[0,50],nexp)");
w.factory("Gaussian:constraint(b0[0,10],b,sigmab[1])");
w.factory("PROD:model(pdf,constraint)");
w.var("b0")->setVal(b);
w.var("b0")->setConstant(true); // needed for being treated as global observables
w.var("sigmab")->setVal(sigmab*b);
ModelConfig mc("ModelConfig",&w);
mc.SetPdf(*w.pdf("model"));
mc.SetParametersOfInterest(*w.var("s"));
mc.SetObservables(*w.var("nobs"));
mc.SetNuisanceParameters(*w.var("b"));
// these are needed for the hypothesis tests
mc.SetSnapshot(*w.var("s"));
mc.SetGlobalObservables(*w.var("b0"));
mc.Print();
// import model in the workspace
w.import(mc);
// make data set with the namber of observed events
RooDataSet data("data","", *w.var("nobs"));
// w.var("nobs")->setVal(3);
w.var("nobs")->setVal(0); // 3 changed to 0.......................
data.add(*w.var("nobs") );
// import data set in workspace and save it in a file
w.import(data);
w.Print();
TString fileName = "CountingModel.root";
// write workspace in the file (recreate file if already existing)
w.writeToFile(fileName, true);
}
and the command is:
StandardHypoTestInvDemo("CountingModel.root", "w", "ModelConfig", "", "data", 0, 2, false, 10, 0, 15, 1000,true )
as you can see, the results is 1.51 but the FC gives the 2.33.
