HistFactory with control- and signal-region

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.

Hi,
The code seems fine, but I would need you to post also the histogram and the full code attached as a file to be able to investigate the issue.
Thanks
Lorenzo

Alright thank you verry much!
I think the problem lies then in my expectation.

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.