Dear stats experts,
I have the following situation at hand:
There is a controlregion (CR) which is completely devoid of any signal which fills a 6 bin counting histogram hobsCR. This histogram is supposed to be transfered to the signalregion (SR) by a transferfactor (of order 0.1) with a systematic error. Now I want to derive an upper limit on a scalefactor mu, which scales the MC signal expectation in the SR given the observed events in the histogram hobsSR.
To generate the model file I use the following:
´ //Create histogram templates for the background
TH1 *h1_bSR = (TH1 *)hobsCR->Clone(“h1_bSR”);
TH1 *h1_bCR = (TH1 *)hobsCR->Clone(“h1_bCR”);HistFactory::Measurement meas(“HistModel”, “HistModel”);
meas.SetPOI(“mu”);meas.SetLumi(1.0);
meas.SetLumiRelErr(0.001);
meas.AddConstantParam(“Lumi”);printf(“make channel SR\n”);
HistFactory::Channel channelSR(“Signalregion”);
channelSR.SetData(hobsSR);//Create signal sample
HistFactory::Sample signal(“signal”);
signal.AddNormFactor(“mu”, 1, 0, 2);
signal.SetHisto(hSignal);
channelSR.AddSample(signal);//Backgroundestimate in SR is CR scaled by the transferfactor
HistFactory::Sample backgSR(“bA”);
backgSR.AddNormFactor(“TransferFac”, transferFac, transferFac - 10 * transferFacErr, transferFac + 10 * transferFacErr);
meas.AddConstantParam(“TransferFac”);
backgSR.AddOverallSys(“transFacConstraint”, 1 - transferFacErr / transferFac, 1 + transferFacErr / transferFac);
backgSR.SetHisto(h1_bSR);
channelSR.AddSample(backgSR);//Create channel for CR
HistFactory::Channel channelCR(“Bregion”);
channelCR.SetData(hobsCR);//Sample in CR
HistFactory::Sample backgB(“bB”);
backgB.SetHisto(h1_bCR);
backgB.ActivateStatError();
channelCR.AddSample(backgB);meas.AddChannel(channelSR);
meas.AddChannel(channelCR);meas.Print();
// make the model
RooWorkspace *w = HistFactory::MakeModelAndMeasurementFast(meas);RooStats::ModelConfig *sbModel = (RooStats::ModelConfig *)w->obj(“ModelConfig”);
RooRealVar *poi = (RooRealVar *)sbModel->GetParametersOfInterest()->first();
ModelConfig *bModel = (ModelConfig *)sbModel->Clone();
bModel->SetName(TString(sbModel->GetName()) + TString(“_with_poi_0”));
double oldval = poi->getVal();
poi->setVal(0);
bModel->SetSnapshot(*poi);
poi->setVal(oldval);
w->import(*bModel);
w->import(*sbModel);w->writeToFile(“HistFactory.root”);
w->Print();`
Using this model I used StandardHypoTestInvDemo with the asymptotic calculator and the one sided profile likelyhood to derive the upper bound on the signal strength mu.
Would that be the correct way to build the model? I doubt my results because the 95% upper limit on mu seems way to strict. The last 3 bins in hobsSR and hobsCR contain fairly few counts (<10) so their poisson errors should allow a lot of signal in this region.
Sorry for the basic question. I have no fluent roostats user in my group.