using namespace RooFit;

void FinalMassFitToysSimple(){

  TChain *chain = new TChain("modelData","");
  chain->Add("ToyDataset_2016.root");
  chain->Add("ToyDataset_2017.root");
  chain->Add("ToyDataset_2018.root");

  RooRealVar x("m","m_{inv} (MeV/c^{2})",3400,3800);
  RooDataSet data("data","m_{inv} (MeV/c^{2})",x,Import(*chain));

  RooRealVar mu("mu","mu",3400,3800);
  mu.setVal(3576.);
  mu.setConstant(kTRUE);
  RooRealVar sig("sig","sig",6.5,6.0,10.0);
  sig.setVal(9.);
  sig.setConstant(kTRUE);
  RooRealVar lambda("lambda", "lambda",  0.001, -0.05, 0.05);

  RooPlot* xframe = x.frame(Title(""));
  data.plotOn(xframe,Name("data"),Binning(40));

  RooGaussian signal("signal", "signal", x, mu, sig);
  RooExponential bkg("bkg", "exponential PDF", x, lambda);
  RooRealVar sigFrac("sigFrac", "Signal Fraction", 0.1, 0, 1);
  RooAddPdf model("model", "Signal plus background PDF",  signal, bkg, sigFrac);

  RooFitResult* resultBkg = bkg.fitTo(data,Save());
  RooFitResult* resultSigBkg = model.fitTo(data,Save());

  model.plotOn(xframe,Name("bkg"),Components(bkg),LineWidth(2),LineColor(kBlue), LineStyle(ELineStyle::kDashed));
  model.plotOn(xframe,Name("signal"),Components(signal),LineWidth(2),LineColor(kRed));
  model.plotOn(xframe,Name("model"));
  xframe->SetTitle("");
  xframe->Draw();

  Double_t qValue;
  Double_t pValue;
  qValue = 2*(resultBkg->minNll() - resultSigBkg->minNll()) > 0.0001 ? 2*(resultBkg->minNll() - resultSigBkg->minNll()) : 0;
  pValue = 0.5*(1. - TMath::Erf(sqrt(qValue/2.)));

  printf("lambda is %f ± %f\n",lambda.getVal(),lambda.getError());
  printf("mean is %f\n",mu.getVal());
  printf("p-value is %.12f\n",pValue);
  printf("significance is %f\n",RooStats::PValueToSignificance(pValue));
  printf("----------------------------------------\n");
  printf("combined fit bkg model minNll is %f\n",resultBkg->minNll());
  printf("combined fit sig+bkg model minNll is %f\n", resultSigBkg->minNll());
}