using namespace RooFit;
using namespace RooStats;

void Priors(double nobs = 5) {
  // test of different priors for Poisson counting experiment

  // number of expected events
  RooRealVar s("s","s",1,0,nobs*3+5);
  // number of observed events
  RooRealVar n("n","n",nobs,0,nobs*2);

  // the likelihood
  RooPoisson pois("pois","Poisson Likelihood",n,s);

  // Priors
  // uniform prior
  RooUniform unifPrior("unifPrior","Uniform prior",s);
  // Jeffreys prior
  RooGenericPdf jefPrior("jefPrior","Jeffreys prior","1.0/sqrt(s)",RooArgList(s));

  // Posteriors
  // Using RooProdPdf
  // uniform prior case
  RooProdPdf unifPrior_posterior("unifPrior_posterior","Posterior for uniform prior",pois,unifPrior);
  // Jeffreys prior case
  RooProdPdf jefPrior_posterior("jefPrior_posterior","Posterior for Jeffreys prior",pois,jefPrior);

  // let's compare the two
  // plot for s
  RooPlot* sframe = s.frame();
  // let's plot Jeffreys case first
  jefPrior_posterior.plotOn(sframe,LineColor(kRed),RooFit::Name("jefPrior_posterior"));
  // then uniform prior one
  unifPrior_posterior.plotOn(sframe,RooFit::Name("unifPrior_posterior"));
  // and then let's draw the result
  sframe->SetTitle(";Signal yield s;Posterior value");
  sframe->Draw();


  // cosmetics
  gPad->SetTicks(1,1);
  // legend
  TString header;
  header.Form("n_{obs} = %d",int(nobs));
  TLegend *leg = new TLegend(0.6,0.6,0.8,0.85,header.Data());
  leg->Draw();
  leg->AddEntry(sframe->findObject("unifPrior_posterior"),"Uniform Prior","l");
  leg->AddEntry(sframe->findObject("jefPrior_posterior"),"Jeffrey Prior","l");

  // let's save this nice result
  TString canvas_name;
  canvas_name.Form("Poisson_posterior_for_%d_nobs.pdf",int(nobs));
  gPad->SaveAs(canvas_name.Data());

  // let's save everything in a workspace
  RooWorkspace w("w");
  // importing pdfs
  w.import(unifPrior);
  w.import(jefPrior);
  w.import(pois);
  w.writeToFile("myPriors.root");

}