{
    using namespace RooFit;
    using namespace RooStats;
    
    RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
    
    TCanvas *c1 = new TCanvas("c1","",800,600);
    c1->SetMargin(0.05,0.05,0.22,0.05);
    c1->Divide(4,4,0.,0.);
    
    printf("\nObserved    0    1    2    3\n");
    for (int b_set=0; b_set<4; b_set++) {
        printf("bkg = %d  ",b_set);
        for (int n_set=0; n_set<4; n_set++) {

            RooRealVar x("x","dummy obs",0.,1.);
            RooUniform pdf_x("pdf_x","dummy pdf of x",x);
            
            RooRealVar s("s","# of signal",1E-5,15.);
            RooRealVar b("b","# of background",b_set);
            RooAddPdf model("model","total PDF",RooArgList(pdf_x,pdf_x),RooArgList(s,b));
            RooUniform prior_s("prior_s","prior for signal",s);
            
            RooWorkspace wspace("wspace");
            ModelConfig cfg(&wspace);
            cfg.SetPdf(model);
            cfg.SetParametersOfInterest(s);
            cfg.SetPriorPdf(prior_s);
            
            RooRealVar nevt("nevt","# of observed event",n_set);
            RooDataSet data("data","data",RooArgSet(x,nevt),WeightVar("nevt"));
            data.add(RooArgSet(x),n_set);
        
            BayesianCalculator bcalc(data,cfg);
            bcalc.SetLeftSideTailFraction(0.);
            bcalc.SetConfidenceLevel(0.9);
            SimpleInterval* interval = bcalc.GetInterval();
            printf("%.2f ",interval->UpperLimit());
            
            c1->cd(b_set*4+n_set+1);
            RooPlot *frame = bcalc.GetPosteriorPlot();
            frame->SetTitle("");
            frame->GetXaxis()->SetTitleSize(0.12);
            frame->GetXaxis()->SetLabelSize(0.10);
            frame->GetYaxis()->SetLabelSize(0.);
            frame->GetYaxis()->SetTitle("");
            frame->Draw();
        }
        printf("\n");
    }
}