#include "StandardHypoTestInvDemo.C"

void CountingModel(){
    Double_t 
    nobs   = 0.602,            // number of observed events
    s      = 0.000,            // number of signal events
    b      = 0.602,            // number of background events
    sigmab = 0.053/0.602;      // relative uncertainty in b

  RooWorkspace w("w");
   
// make Poisson model * Gaussian constraint
w.factory("sum:nexp(s[0,0,10],b[1,0,10])");
// Poisson of (n | s+b)
w.factory("Poisson:pdf(nobs[0,30],nexp)");
w.factory("Gaussian:constraint(b0[0,10],b,sigmab[0.5])");
w.factory("PROD:model(pdf,constraint)");

w.var("b")->setVal(b);
w.var("nobs")->setVal(nobs); 
w.var("s")->setVal(s);


w.var("b0")->setVal(b);
w.var("b0")->setConstant(true); // needed for being treated as global observables
w.var("sigmab")->setVal(sigmab*b);  


ModelConfig mc("ModelConfig",&w);
mc.SetPdf(*w.pdf("model"));
mc.SetParametersOfInterest(*w.var("s"));
mc.SetObservables(*w.var("nobs"));
mc.SetNuisanceParameters(*w.var("b"));

// these are needed for the hypothesis tests
mc.SetSnapshot(*w.var("s"));
mc.SetGlobalObservables(*w.var("b0"));

mc.Print();
// import model in the workspace 
w.import(mc);

// make data set with the namber of observed events
RooDataSet dtst("dtst","", *w.var("nobs"));
//  w.var("nobs")->setVal(3);
 

dtst.add(*w.var("nobs") );
// import data set in workspace and save it in a file
w.import(dtst);

w.Print();

w.writeToFile("CountingModel.root", true);

  StandardHypoTestInvDemo(
                          "CountingModel.root",   // file name
                          "w",                    // workspace name
                          "ModelConfig",          // S+B modelconfig name
                          "",                     // B model name
                          "dtst",                 // data set name
                          0,                      // calculator type (0 = Freq)
                          2,                      // test statistic type (2 = two-sided Profile Likelihood)
                          false,                   // use CLs
                          25,                     // number of points
                          0,                      // xmin
                          4,                      // xmax
                          1e3                     // number of toys
                          ); 
}