#include <RooWorkspace.h>
#include <RooRealVar.h>
#include <RooAbsPdf.h>
#include <RooMinuit.h>
#include <RooDataSet.h>
#include <RooPlot.h>

#include <TCanvas.h>
#include <TPad.h>
#include <TAxis.h>

void likelihood () {
  // Construct a model
  RooWorkspace w("w", kTRUE);

  w.factory("Gaussian::g1(x[-20, 20], mean[-10, 10], sigma_g1[3])");
  w.factory("Gaussian::g2(x, mean, sigma_g2[4, 3, 6])");
  w.factory("SUM::model(frac[0.5, 0, 1] * g1, g2)");

  // Generate 1000 events
  RooDataSet *data = w.pdf("model")->generate(*w.var("x"), 1000);

  // Create a likelihood function
  RooAbsReal *nll = w.pdf("model")->createNLL(*data, RooFit::NumCPU(4));

  // Minimize the likelihood
  RooMinuit(*nll).migrad();

  // Plot the likelihood scan in parameter frac
  RooPlot *frame1 = w.var("frac")->frame(RooFit::Bins(10), RooFit::Range(0.01, 0.95));
  nll->plotOn(frame1, RooFit::ShiftToZero());

  // Plot the profile likelihood in frac
  RooAbsReal *pll_frac = nll->createProfile(*w.var("frac"));
  pll_frac->plotOn(frame1, RooFit::LineColor(kRed));

  // Plot the likelihood scan in parameter sigma_g2
  RooPlot *frame2 = w.var("sigma_g2")->frame(RooFit::Bins(10), RooFit::Range(3.3, 5.0));
  nll->plotOn(frame2, RooFit::ShiftToZero());

  // Plot the profile likelihood in sigma_g2
  RooAbsReal *pll_g2 = nll->createProfile(*w.var("sigma_g2"));
  pll_g2->plotOn(frame2, RooFit::LineColor(kRed));

  // Draw all frames on a canvas
  TCanvas *c = new TCanvas("likelihood", "likelihood", 800, 400);
  c->Divide(2);
  c->cd(1);
  gPad->SetLeftMargin(0.15);
  frame1->GetYaxis()->SetTitleOffset(1.6);
  frame1->Draw();
  c->cd(2);
  gPad->SetLeftMargin(0.15);
  frame2->GetYaxis()->SetTitleOffset(1.6);
  frame2->Draw();
} // likelihood