void test() {
  using namespace RooFit;
  using namespace RooStats;

  double CL   = 0.95;
  RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR);

  RooWorkspace* wspace = new RooWorkspace("wspace");
  wspace->factory("NBobs[7861.95]");
  wspace->factory("NTdata[709.00]");
  wspace->factory("NTbck[85.61,0,2000]");
  wspace->var("NTbck")->setVal(85.61);  
  wspace->var("NTbck")->setConstant();
  
  RooRealVar *BR = wspace->factory("BR[0,-0.000001,1.]");

  wspace->factory("Xsys[1,-10,10]");
  wspace->var("Xsys")->setVal(0.);
  //wspace->var("Xsys")->setConstant();

  wspace->factory("CS_tt[157.5,0,300]"); 
  wspace->var("CS_tt")->setVal(157.5);
  wspace->var("CS_tt")->setConstant(); 

  wspace->factory("CS_tW[15.74,0,40]");
  wspace->var("CS_tW")->setVal(15.74);
  wspace->var("CS_tW")->setConstant();

  wspace->factory("expr::effB_SM_SM('exp(log(0.00904)+Xsys*log(1+10.51))',Xsys)");
  wspace->factory("expr::effB_BNV_SM('exp(log(0.06045)+Xsys*log(1+10.51))',Xsys)");
  wspace->factory("expr::effB_BNV_BNV('exp(log(0.03551)+Xsys*log(1+10.51))',Xsys)");
  wspace->factory("expr::effB_SM('exp(log(0.00262)+Xsys*log(1+10.51))',Xsys)");
  wspace->factory("expr::effB_BNV('exp(log(0.01824)+Xsys*log(1+10.51))',Xsys)");

  wspace->factory("expr::effT_SM_SM('exp(log(0.00065)+Xsys*log(1+10.51))',Xsys)");
  wspace->factory("expr::effT_BNV_SM('exp(log(0.01756)+Xsys*log(1+10.51))',Xsys)");
  wspace->factory("expr::effT_BNV_BNV('exp(log(0.00912)+Xsys*log(1+10.51))',Xsys)");
  wspace->factory("expr::effT_SM('exp(log(0.00016)+Xsys*log(1+10.51))',Xsys)");
  wspace->factory("expr::effT_BNV('exp(log(0.00552)+Xsys*log(1+10.51))',Xsys)");

  wspace->factory("expr::effB_TT('2*BR*(1-BR)*effB_BNV_SM+(1-BR)*(1-BR)*effB_SM_SM+BR*BR*effB_BNV_BNV',BR,effB_BNV_SM,effB_SM_SM,effB_BNV_BNV)");
  wspace->factory("expr::effT_TT('2*BR*(1-BR)*effT_BNV_SM+(1-BR)*(1-BR)*effT_SM_SM+BR*BR*effT_BNV_BNV',BR,effT_BNV_SM,effT_SM_SM,effT_BNV_BNV)");

  wspace->factory("expr::effB_TW('(1-BR)*effB_SM+BR*effB_BNV',BR,effB_SM,effB_BNV)");
  wspace->factory("expr::effT_TW('(1-BR)*effT_SM+BR*effT_BNV',BR,effT_SM,effT_BNV)");

  wspace->factory("expr::A('1 / (1 + (CS_tW/CS_tt)*(effB_TW/effB_TT))',CS_tW,CS_tt,effB_TW,effB_TT)");
  wspace->factory("expr::B('1 / (1 + (CS_tt/CS_tW)*(effB_TT/effB_TW))',CS_tt,CS_tW,effB_TT,effB_TW)");
  
  wspace->factory("expr::NTexp('NBobs*(A*(effT_TT/effB_TT)+B*(effT_TW/effB_TW))+NTbck',NBobs,A,effT_TT,effB_TT,B,effT_TW,effB_TW,NTbck)");

  /*
  wspace->factory("expr::NTexpA('NBobs*((1/(1+(CS_tW/CS_tt)*(((1-BR)*(exp(log(0.00262)+Xsys*log(1+10.51)))+BR*exp(log(0.01824)+Xsys*log(1+10.51)))/(2*BR*(1-BR)*(exp(log(0.06045)+Xsys*log(1+10.51)))+(1-BR)*(1-BR)*(exp(log(0.00904)+Xsys*log(1+10.51)))+BR*BR*(exp(log(0.03551)+Xsys*log(1+10.51)))))))*((2*BR*(1-BR)*(exp(log(0.01756)+Xsys*log(1+10.51)))+(1-BR)*(1-BR)*exp(log(0.00065)+Xsys*log(1+10.51))+BR*BR*(exp(log(0.00912)+Xsys*log(1+10.51))))/(2*BR*(1-BR)*(exp(log(0.06045)+Xsys*log(1+10.51)))+(1-BR)*(1-BR)*(exp(log(0.00904)+Xsys*log(1+10.51)))+BR*BR*(exp(log(0.03551)+Xsys*log(1+10.51))))))',NBobs,CS_tW,CS_tt,BR,Xsys)");

  wspace->factory("expr::NTexpB('NBobs*((1/(1+(CS_tt/CS_tW)*((2*BR*(1-BR)*(exp(log(0.06045)+Xsys*log(1+10.51)))+(1-BR)*(1-BR)*(exp(log(0.00904)+Xsys*log(1+10.51)))+BR*BR*(exp(log(0.03551)+Xsys*log(1+10.51))))/((1-BR)*(exp(log(0.00262)+Xsys*log(1+10.51)))+BR*(exp(log(0.01824)+Xsys*log(1+10.51)))))))*(((1-BR)*(exp(log(0.00016)+Xsys*log(1+10.51)))+BR*(exp(log(0.00552)+Xsys*log(1+10.51))))/((1-BR)*(exp(log(0.00262)+Xsys*log(1+10.51)))+BR*(exp(log(0.01824)+Xsys*log(1+10.51))))))+NTbck',NBobs,CS_tW,CS_tt,BR,Xsys,NTbck)"); 
  */

  RooRealVar *dataVal = wspace->var("NTdata");
  dataVal->setVal(709.00);
  RooDataSet *NTdat = new RooDataSet("data","data",*dataVal);
  NTdat->add(*dataVal);
  wspace->import(*data);

  RooAbsPdf* model = (RooAbsPdf*) wspace->factory("Poisson::model_(NTdata,NTexp)"); //sum(NTexpA,NTexpB))");
  wspace->defineSet("obs","NTdata");
  wspace->defineSet("parOfInterest", "BR"); 
  wspace->defineSet("nuisPar", "Xsys");
  
  wspace->factory("Gaussian::sys(Xsys,0,1)");
  //wspace->factory("Gaussian::sys(NTbck,85,40)");
  wspace->factory("PROD::model(model_,sys)");

  ModelConfig* modelConfig = new ModelConfig("SimpleCount");  
  modelConfig->SetWorkspace(*wspace);
  modelConfig->SetPdf(*wspace->pdf("model"));
  modelConfig->SetParametersOfInterest(*wspace->set("parOfInterest")); 
  modelConfig->SetNuisanceParameters(*wspace->set("nuisPar"));
  wspace->import(*modelConfig); 
  
  ProfileLikelihoodCalculator plc(*data, *modelConfig); 
  RooArgSet* nullParams = (RooArgSet*)(*wspace->set("parOfInterest")).snapshot(); 
  nullParams->setRealValue("BR",0); 
  plc.SetNullParameters(*nullParams); 
  HypoTestResult* htr = plc.GetHypoTest();
  double pValue = htr->NullPValue();
  double significance = TMath::NormQuantile(1 - pValue); 

  plc.SetConfidenceLevel(0.95);
  LikelihoodInterval* plInt = plc.GetInterval();
  LikelihoodIntervalPlot* intPlot = new LikelihoodIntervalPlot(plInt);
  double Lo_pl = plInt->LowerLimit(*BR);
  double Up_pl = plInt->UpperLimit(*BR);
 
  cout << "Significance = " << significance << endl;                
  cout << "Lower limit  = " << Lo_pl << endl;
  cout << "Upper limit  = " << Up_pl <<endl;
  gApplication->Terminate();  
}