Thanks for the reply. I have an additional question about it. When hypotestcalculator
use CLs method, it should calculate profile likelihood ratio. In that likelihood, is Poisson term included?
Actually, I wanted to make unbinned analysis. Here is core part of my code:
(... read root file and do something ...)
RooWorkspace* w = new RooWorkspace("w", "workspace");
// obs
w->factory("Rarity[0,1]");
// pdfs
w->factory("EXPR::bkg('(1-Rarity)*exp(Rarity*a+Rarity*Rarity*b)',Rarity,a[-1.883],b[1.09])");
w->factory("Polynomial::sig(Rarity)");
// systematics
w->factory("Gaussian::add_sub(0,alpha[-10.0,10.0],1.39)");
w->factory("Gaussian::mul_sub(0,beta[-0.5,0.5],0.0555)");
w->factory("sig_eff[0.0014]");
w->factory("NBB[330000000]");
w->factory("sigBR[0.00008, 0.0, 0.0005]");
w->factory("prod::sig_yield(2,NBB,sig_eff,sigBR)");
w->factory("expr::sig_yield_with_unc('sig_yield*(1+beta)+alpha', sig_yield, alpha, beta)");
w->factory("bkg_yield[1900, 1700, 2100]");
w->factory("SUM::model(sig_yield*sig, bkg_yield*bkg)");
w->factory("PROD::model_with_unc(model,add_sub,mul_sub)");
RooStats::ModelConfig mc("ModelConfig", w);
mc.SetPdf(*w->pdf("model_with_unc"));
mc.SetParametersOfInterest(*w->var("sigBR"));
mc.SetNuisanceParameters(RooArgSet(*w->var("alpha"), *w->var("beta"), *w->var("bkg_yield")));
w->import(*d_Rarity, RooFit::Rename("observed_data"));
w->import(mc);
w->Print();
w->writeToFile("CLs_test.root");
, where alpha
and beta
are additive and multiplicative systematic uncertainties, respectively. Signal pdf is flat function (0th order polynomial). Background pdf is slightly complicated function ((1-Rarity)*exp(Rarity*a+Rarity*Rarity*b
). When I use hypotestcalculator, I did
(... read workspace...)
RooStats::ModelConfig* sbModel = (RooStats::ModelConfig*)w->obj("ModelConfig");
RooStats::ModelConfig* bModel = (RooStats::ModelConfig*)sbModel->Clone("BonlyModel");
RooRealVar* poi = (RooRealVar*)bModel->GetParametersOfInterest()->first();
poi->setVal(0);
bModel->SetSnapshot(*poi);
RooStats::FrequentistCalculator hybCalc(*data, *bModel, *sbModel);
RooStats::ProfileLikelihoodTestStat* plr = new RooStats::ProfileLikelihoodTestStat(*sbModel->GetPdf());
plr->SetOneSided(true);
RooStats::ToyMCSampler* toymcs = (RooStats::ToyMCSampler*)hybCalc.GetTestStatSampler();
toymcs->SetTestStatistic(plr);
hybCalc.SetToys(1000, 1000);
(... do something...)
Then, I can obtain upper limit from CLs method. Does profile likelihood, which is used for CLs method, include Poisson term? ((e^(-\mu))*\mu^N/N!
)
Thank you!