void FillNonzero(TH1* h) {
  for(int i = 1; i <= h->GetNbinsX(); ++i) {
    if(h->GetBinContent(i) == 0) h->SetBinContent(i, 1e-10);
  }
}

void pdf2(const int kN1, const int kN2, const double kP1) {
  const double kXmin = -3;
  const double kXmax =  10;
  const int kNbins = 130;
  const double kDx = 3;

  gRandom->SetSeed(1000);

  TH1D* h1 = new TH1D("h1", "h1", kNbins, kXmin, kXmax);
  h1->FillRandom("gaus", kN1);
  FillNonzero(h1);

  TH1D* h2 = new TH1D("h2", "h2", kNbins, kXmin, kXmax);
  h2->FillRandom("gaus", kN2 * (1 - kP1));

  TF1 fP1("fP1", Form("exp(-(x - %f)**2/2.)", kDx), kXmin, kXmax);
  h2->FillRandom("fP1", kN2 * kP1);

  TF1 fP2("fP2", Form("exp(-(x - %f)**2/2.)", kDx * 2), kXmin, kXmax);
  h2->FillRandom("fP2", kN2 * kP1 * 0.5);

  TF1 fP3("fP3", Form("exp(-(x - %f)**2/2.)", kDx * 3), kXmin, kXmax);
  h2->FillRandom("fP3", kN2 * kP1 * 0.2);

  RooRealVar x("x", "x", kXmin, kXmax);
  RooRealVar p0("p0", "p0", 0.5 * kN2, 0., 1.5 * kN2);
  RooRealVar p1("p1", "p1", 0.5 * kN2, 0., 1.5 * kN2);
  RooRealVar p2("p2", "p2", 0.5 * kN2, 0., 1.5 * kN2);

  RooDataHist dh1("dh1", "dh1", x, h1);
  RooDataHist dh2("dh2", "dh2", x, h2);

  const int kIntOrder = 0; // 1 for interpolation

  RooHistPdf pdf0("pdf0", "pdf0", x, dh1, kIntOrder);

  RooRealVar dx("dx", "dx", 2.5, 2., 4.);
  RooFormulaVar dx1("dx1", "x - dx", RooArgSet(x, dx));
  RooHistPdf pdf1("pdf1", "pdf1", dx1, x, dh1, kIntOrder);

  RooFormulaVar dx2("dx2", "x - 2*dx", RooArgSet(x, dx));
  RooHistPdf pdf2("pdf2", "pdf2", dx2, x, dh1, kIntOrder);

  RooAddPdf model("model", "model", RooArgList(pdf0, pdf1, pdf2), RooArgList(p0, p1, p2));

  RooAbsReal* nll = model.createNLL(dh2, RooFit::NumCPU(8));
  RooMinuit(*nll).migrad();

  TCanvas* c1 = new TCanvas("c1", "c1");
  c1->Divide(2, 2);
  RooPlot* frameNLL = dx.frame(RooFit::Bins(1000), RooFit::Range(dx.getValV() * 0.5, dx.getValV() * 1.5));
  nll->plotOn(frameNLL, RooFit::ShiftToZero());
  c1->cd(1);
  frameNLL->Draw();
  frameNLL->GetXaxis()->SetRangeUser(dx.getValV() * 0.8, dx.getValV() * 1.2);
  frameNLL->GetYaxis()->SetRangeUser(0, 200);

  frameNLL = p0.frame(RooFit::Bins(1000), RooFit::Range(p0.getValV() * 0.8, p0.getValV() * 1.2));
  nll->plotOn(frameNLL, RooFit::ShiftToZero());
  c1->cd(2);
  frameNLL->Draw();
  frameNLL->GetXaxis()->SetRangeUser(p0.getValV() * 0.95, p0.getValV() * 1.05);
  frameNLL->GetYaxis()->SetRangeUser(0, 5);

  frameNLL = p1.frame(RooFit::Bins(1000), RooFit::Range(p1.getValV() * 0.5, p1.getValV() * 1.5));
  nll->plotOn(frameNLL, RooFit::ShiftToZero());
  c1->cd(3);
  frameNLL->Draw();
  frameNLL->GetXaxis()->SetRangeUser(p1.getValV() * 0.9, p1.getValV() * 1.1);
  frameNLL->GetYaxis()->SetRangeUser(0, 5);

  frameNLL = p2.frame(RooFit::Bins(1000), RooFit::Range(p2.getValV() * 0.5, p2.getValV() * 1.5));
  nll->plotOn(frameNLL, RooFit::ShiftToZero());
  c1->cd(4);
  frameNLL->Draw();
  frameNLL->GetXaxis()->SetRangeUser(p2.getValV() * 0.9, p2.getValV() * 1.1);
  frameNLL->GetYaxis()->SetRangeUser(0, 5);

  TCanvas* c2 = new TCanvas("c2", "c2");
  RooPlot* xframe = x.frame(RooFit::Title("Gaussian Fit Test"));
  dh2.plotOn(xframe);
  model.plotOn(xframe);
  model.plotOn(xframe, RooFit::Components(pdf0), RooFit::LineStyle(ELineStyle::kDashed));
  model.plotOn(xframe, RooFit::Components(pdf1), RooFit::LineStyle(ELineStyle::kDashed));
  model.plotOn(xframe, RooFit::Components(pdf2), RooFit::LineStyle(ELineStyle::kDashed));

  xframe->Draw();
}