/*
    Interactive exercise:
    Can vary the uncertainty on the background?
    Can you reproduce the results? Why?
    What changes once binning? Why?

    # constant Bkg
    Limits: [39.6066 - 80.4865]
    Significance (SL2): 3.07742

    # 0.1% Sys on Bkg
    Limits: [39.6068 - 80.4931]
    Significance (SL2): 3.07703

    # 30% Sys on Bkg
    Limits: [40.5089 - 90.4514]
    Significance (SL2): 2.67559
    */

    {
    using namespace RooFit; using namespace RooStats;

    // Model built with Low Level Factory
    RooWorkspace w("CMGWs","CMG workspace");
    w.factory("Gaussian::signal(imass[110,120],sigmean[115,105,120],sigwidth[1])");
    w.factory("Exponential::background(imass,exppar[-.1,-1,0])");
    w.factory("SUM::modelbase(nsig[50,0,120]*signal,nbkg[1000,0,10000]*background)");  
    w.factory("Gaussian::nbkgConstraint(nbkg,1000,1)"); // 0.1% background uncertainty
    w.factory("PROD::model(modelbase,nbkgConstraint)"); // implement a systematic

    // Get the interesting objects
    RooAbsPdf* model = w.pdf("model");
    RooRealVar* obs = w.var("imass"); obs->setBins(30);
    RooRealVar* nsig = w.var("nsig");RooRealVar* nbkg = w.var("nbkg");

    // Create an example dataset (binned is fast..)
    RooDataHist* data = model->generateBinned(*obs);

    // Adapt to fluctuations
    model->fitTo(*data);

    // Model Config
    ModelConfig modelConfig(new RooWorkspace());
    modelConfig.SetPdf(*model);
    modelConfig.SetParametersOfInterest(RooArgSet(*nsig));
      
    // Let's make a plot and some maquillage
    TCanvas* dataCanvas = new TCanvas("res","The Result",800,400);
    dataCanvas->Divide(2,1); dataCanvas->cd(1);

    RooPlot* imassframe = obs->frame();
    data->plotOn(imassframe); model->plotOn(imassframe);
    imassframe->Draw();

    // Let's use a Calculator based on the Profile Likelihood Ratio
    double CL = .68; 
    double size = 1-CL; // size of the test

    ProfileLikelihoodCalculator plc(*data,modelConfig,size); 
    LikelihoodInterval* lrint = plc.GetInterval(); 
    LikelihoodIntervalPlot plotInt(lrint);

    dataCanvas->cd(2); plotInt.Draw();

    // Calculate the significance: SL2 estimator
    nsig->setVal(0);
    plc.SetNullParameters(*nsig);
    HypoTestResult* hypotest = plc.GetHypoTest();
    double SL2 = hypotest->Significance();

    cout << "Limits: [" << lrint.LowerLimit(*nsig) << " - " 
                        << lrint.UpperLimit(*nsig) << "]" << endl;
    cout << "Significance (SL2): " << SL2 << endl;

    }