using namespace RooFit;

void buildModel(RooWorkspace *wspace)
{
    TFile *fin = new TFile("example_06.root");
    TNtupleD* nt = (TNtupleD *)fin->Get("nt");
    
    RooRealVar mass("mass","mass obs",0.,2.);
    
    RooRealVar mu("mu","signal mean",1.0);
    RooRealVar sigma("sigma","signal width",0.05);
    RooGaussian gaus("gaus","signal PDF",mass,mu,sigma);
    
    RooRealVar slope("slope","background slope",-0.3,-5.,5.);
    RooPolynomial linear("linear","background PDF",mass,RooArgSet(slope));
    
    RooRealVar ns("ns","ns",10,0.,1000.);
    RooRealVar nb("nb","nb",90,0.,1000.);
    RooAddPdf model("model","full PDF",RooArgList(gaus,linear),RooArgList(ns,nb));
    
    RooDataSet data("data","data",nt,RooArgSet(mass));
    
    wspace->import(model);
    wspace->import(data);

    delete fin;
}

void example_07()
{
    RooWorkspace *wspace = new RooWorkspace("wspace");
    buildModel(wspace);
    
    RooFitResult *res0 = wspace->pdf("model")->fitTo(*wspace->data("data"),Save(true),Minos(true));
    
    TH1D *scan_2nll = new TH1D("scan_2nll","",101,-0.1,20.1);
    TH1D *scan_beta = new TH1D("scan_beta","",101,-0.1,20.1);
    
    for (int i=1; i<=scan_2nll->GetNbinsX(); i++) {
        wspace->var("ns")->setVal(scan_2nll->GetBinCenter(i));
        wspace->var("ns")->setConstant(true);
        RooFitResult *res1 = wspace->pdf("model")->fitTo(*wspace->data("data"),Save(true));
        double d2NLL = (res1->minNll()-res0->minNll())*2.;
        
        scan_2nll->SetBinContent(i,d2NLL);
        scan_beta->SetBinContent(i,1.-TMath::Prob(d2NLL,1));
        
        delete res1;
    }
    
    RooRealVar* best_ns = (RooRealVar*)res0->floatParsFinal().find("ns");
    
    TCanvas *c1 = new TCanvas("c1","",800,400);
    c1->Divide(2);
    
    TLine lin;
    lin.SetLineColor(kGray+2);
    lin.SetLineStyle(kDashed);
    
    TBox box;
    box.SetFillColor(kRed-10);
    
    c1->cd(1);
    c1->GetPad(1)->SetMargin(0.14,0.06,0.13,0.07);
    
    scan_2nll->SetStats(false);
    scan_2nll->SetLineWidth(2);
    scan_2nll->GetXaxis()->SetTitle("Signal Yield");
    scan_2nll->GetYaxis()->SetTitle("-2ln(L/L_{max})");
    scan_2nll->GetXaxis()->SetTitleSize(0.06);
    scan_2nll->GetYaxis()->SetTitleSize(0.06);
    scan_2nll->GetYaxis()->SetTitleOffset(1.0);
    scan_2nll->Draw("axis");
    
    box.DrawBox(best_ns->getVal()+best_ns->getErrorLo(),0.,
                best_ns->getVal()+best_ns->getErrorHi(),scan_2nll->GetMaximum());
    
    for(int n=0; n<=2; n++)
        lin.DrawLine(0.,n*n,20.,n*n);
    
    scan_2nll->Draw("csame");
    
    c1->cd(2);
    c1->GetPad(2)->SetMargin(0.14,0.06,0.13,0.07);
    
    scan_beta->SetStats(false);
    scan_beta->SetLineWidth(2);
    scan_beta->GetXaxis()->SetTitle("Signal Yield");
    scan_beta->GetYaxis()->SetTitle("Confidence Level (#beta)");
    scan_beta->GetXaxis()->SetTitleSize(0.06);
    scan_beta->GetYaxis()->SetTitleSize(0.06);
    scan_beta->GetYaxis()->SetTitleOffset(1.0);
    scan_beta->Draw("axis");
    
    box.DrawBox(best_ns->getVal()+best_ns->getErrorLo(),0.,
                best_ns->getVal()+best_ns->getErrorHi(),scan_beta->GetMaximum());
    
    for(auto prob: {0.,0.6827,0.9545})
        lin.DrawLine(0.,prob,20.,prob);
    
    scan_beta->Draw("csame");
}