Dear RooStats users,
I am quite new to the RooStats software, and I am trying to set a 90% CL upper limit on a signal strength parameter, that enters in an “on-off” problem.
My PDF reads:
P(n,n’;s,b) = Poisson(n,s+b) * Posson (n’,b)
where n and n’ are my observables, s is the POI, and b a nuisance parameter. I am simplifying the problem here, assuming that the ancillary measurement to determine b has the same luminosity / efficiency as the main one.
I have a doubt regarding how to handle n’. I see these two options:
OPTION-1: handle n’ as a “global fixed variable”
int nsb = 30; // number of observed events
int nb = 28; // number of background events
RooWorkspace w("w",true);
// make Poisson model
w.factory("sum:nexp(s[2,0,20],b[0,0,50])");
// Poisson of (nsb_obs | s+b)
w.factory("Poisson:psb(nsb_obs[50,0,70],nexp)");
// Poisson of (nb_obs | b)
w.factory("Poisson:pb(nb_obs[50,0,50],b)");
//The model
w.factory("PROD:model(psb,pb)");
w.var("nb_obs")->setVal(nb);
w.var("nb_obs")->setConstant(true); // needed for being treated as global observables
//set value of observed events
RooRealVar *nsb_obs = w.var("nsb_obs");
nsb_obs->setVal(nsb);
// make data set with the number of observed events
RooDataSet data("data","",*nsb_obs);
data.add(*nsb_obs);
w.import(data);
w.Print();
ModelConfig mc("Smodel",&w);
mc.SetPdf(*w.pdf("model"));
mc.SetParametersOfInterest(*w.var("s"));
mc.SetObservables(*w.var("nsb_obs"));
// these are needed for the hypothesis tests
mc.SetSnapshot(*w.var("s"));
mc.SetGlobalObservables(*w.var("nb_obs"));
OPTION-2: handle n’ as an observable
// make Poisson model
w.factory("sum:nexp(s[2,0,20],b[0,0,50])");
// Poisson of (nsb_obs | s+b)
w.factory("Poisson:psb(nsb_obs[50,0,70],nexp)");
// Poisson of (nb_obs | b)
w.factory("Poisson:pb(nb_obs[50,0,50],b)");
//The model
w.factory("PROD:model(psb,pb)");
//set value of observed events
RooRealVar *nsb_obs = w.var("nsb_obs");
nsb_obs->setVal(nsb);
//set value of observed events in ancillary measurement
RooRealVar *nb_obs = w.var("nb_obs");
nb_obs->setVal(nb);
// make data set with the number of observed events in both measurements
RooDataSet data("data","", RooArgSet(*nsb_obs,*nb_obs ));
data.add(RooArgSet(*nsb_obs,*nb_obs ));
w.import(data);
w.Print();
ModelConfig mc("Smodel",&w);
mc.SetPdf(*w.pdf("model"));
mc.SetParametersOfInterest(*w.var("s"));
mc.SetObservables(RooArgSet(*w.var("nsb_obs"),*w.var("nb_obs")));
mc.SetSnapshot(*w.var("s"));
With both these strategies, I can compute the 90% CL limit, using for example the StandardHypoTestInvDemo.C (I omitted the code to create the B-only model).
I get same results concerning the CL limit (CLs+b -> s=11.8464), but different values for the expected upper limits, using the B (alternate) model.
What is the difference between these two approaches? I think that, in this case, nb_obs should be declared as a “observables”, i.e. as a probabilistic variables?
Thanks,
Bests,
Andrea