using namespace RooFit;
using namespace RooStats;


void jan_2spinAsy_constrain(double nObs1=2, double nObs2=4, double Amax=0.6){
 
  RooWorkspace* ws = new RooWorkspace("jan2Asy");

  // variables:  predicetd  yields
  RooRealVar* n1 = (RooRealVar*) ws->factory("n1[15,0,50]"); 
  RooRealVar* n2 = (RooRealVar*) ws->factory("n2[15,0,50]"); 

  // model parameterization
  RooRealVar*  A= (RooRealVar*)  ws->factory(Form("A[0.12,%.2f,%.2f]",-Amax,Amax));
  RooRealVar* N0= (RooRealVar*)  ws->factory("N0[15,0,50]");

  /* decide what role plays each variable */  
  ws->defineSet("obs","n1,n2");
  ws->defineSet("parOfInterest", "A");
  ws->defineSet("nuisPar", "N0");
  
  // intermediate variables depend on the true params of the model
  ws->factory("expr::mu1('N0*(1+A)',N0,A)");
  ws->factory("expr::mu2('N0*(1-A)',N0,A)"); 
 
  // define pdf's
  RooAbsPdf* model1 = (RooAbsPdf*) ws->factory("Poisson::pdf1(n1,mu1)");
  RooAbsPdf* model2 = (RooAbsPdf*) ws->factory("Poisson::pdf2(n2,mu2)");
  ws->factory("PROD::pdfTotal(pdf1,pdf2)"); 


  // First create empty data set, set data value and put into workspace
  RooDataSet* dataS = new RooDataSet("dataJan", "data 2spin", RooArgSet(*n1,*n2));  // empty
  // add result of experiment: a pair of 2 measurements for + & - spin state
  n1->setVal(nObs1);   n2->setVal(nObs2); 
  dataS->add(RooArgSet(*n1,*n2));
  dataS->Print("v");
  ws->import(*dataS);
  ws->Print();

  // return;

  // Make ModelConfig object and put in workspace
  ModelConfig* modelConfig = new ModelConfig("recoSSA");  
  modelConfig->SetWorkspace(*ws);
  modelConfig->SetPdf("pdfTotal");
  modelConfig->SetParametersOfInterest(*ws->set("parOfInterest"));
  modelConfig->SetNuisanceParameters(*ws->set("nuisPar"));
  ws->import(*modelConfig);

 
  ProfileLikelihoodCalculator plc(*dataS, *modelConfig);
  
  double CL=0.6827;
  plc.SetConfidenceLevel(CL);
  LikelihoodInterval* plInt = plc.GetInterval();
  double Lo = plInt->LowerLimit(*A);
  double Up = plInt->UpperLimit(*A);
  TString outTxt=Form("ProfLike CL=%.4f -->  A limit[%.3f, %.3f]\n",CL,Lo,Up);
  cout<<outTxt<<endl;
  double asyErr=2*sqrt(nObs1*nObs2/pow(nObs1+nObs2,3));// asy=a-b/a+b
  double recoA=(nObs1-nObs2)/(nObs1+nObs2);
  printf(" INPUT:  nObs1=%.1f  nObs2=%.1f -->Gauss method recA=%.3f  1sigma interval[%.3f, %.3f]\nA\n",nObs1,nObs2,recoA,recoA-asyErr,recoA+asyErr); 
  
  // return;
  
  // PLOT  likelihood for 1-sigma

  c=new TCanvas();
  LikelihoodIntervalPlot* intPlot = new LikelihoodIntervalPlot(plInt);
  intPlot->Draw();
  tx=new TText(-0.5,1.6,outTxt); tx->Draw();

}