#include <TFile.h>
#include <TH1.h>
#include <TStopwatch.h>
#include <RooRealVar.h>
#include <RooConstVar.h>
#include <RooDataHist.h>
#include <RooCFunction3Binding.h>
#include <RooCFunction4Binding.h>
#include <RooFormulaVar.h>
#include <RooFFTConvPdf.h>
#include <RooExtendPdf.h>
#include <RooAddPdf.h>
#include <RooArgusBG.h>
#include <RooSimultaneous.h>
#include <RooFitResult.h>

#include <cmath>

using namespace RooFit;
using namespace std;

static TString Mean() { return "m_{#phi}"; }
static TString Sigma() { return "#sigma"; }
static TString Gamma() { return "#Gamma"; }


inline double GetKMass() { return 493.677; }
inline double GetKThreshold() { return 2 * GetKMass(); }

TString GetPolyParName(const int i) { return TString::Format("a%d", i); }


inline
double
qGaussian(const double x, const double mean, const double sigma, const double q)
{
  const double power = 1.0 / (1.0 - q);
  const double alpha = 0.5 * (x - mean) * (x - mean) / sigma / sigma;
  return std::pow(1 + (q - 1) * alpha, power);
}


double
NA49RelativisticBW(const double m, const double mean, const double gamma)
{
  static const double mk = GetKMass();
  static const double th = GetKThreshold();

  if (m <= th)
    return 0;

  const double mk2 = mk * mk;
  const double m2 = m * m;
  const double m02 = mean * mean;

  const double q2 = m2 / 4.0 - mk2;
  const double q02 = m02 / 4.0 - mk2;
  const double q = std::sqrt(q2);
  const double q0 = std::sqrt(q02);

  const double gammaOfM = 2 * gamma * std::pow(q / q0, 3) * q02 / (q2 + q02);

  const double denominator = (m2 - m02) * (m2 - m02) +
      m02 * gammaOfM * gammaOfM;

  return m * gammaOfM / denominator;
}


void
test()
{
  const int nbins = 104;
  const double xmin = 986;
  const double xmax = 1090;
  TH1D h("h", "", nbins, xmin, xmax);
  const double contents[nbins] = {0, 26, 30, 54, 55, 38, 52, 63, 60, 54, 71, 79, 57, 74, 57, 72, 78, 76, 67, 84, 67, 87, 81, 88, 101, 102, 109, 111, 138, 170, 199, 264, 352, 379, 350, 291, 231, 207, 159, 136, 143, 125, 111, 123, 105, 110, 122, 112, 110, 104, 105, 110, 113, 108, 81, 108, 96, 81, 113, 111, 85, 92, 109, 118, 110, 115, 91, 98, 118, 101, 105, 119, 99, 111, 112, 106, 114, 112, 103, 105, 101, 110, 110, 106, 104, 145, 113, 99, 121, 130, 98, 87, 107, 111, 99, 119, 113, 121, 117, 112, 109, 111, 111, 108};
  for (int i = 1; i <= nbins; ++i)
    h.SetBinContent(i, contents[i - 1]);

  RooRealVar x("x", "x", xmin, xmax);
  x.setRange("norm", xmin, xmax);
  x.setRange("cache", xmin, xmax);
  x.setBins(10000, "cache");

  RooDataHist data("data", "data", x, &h);

  RooRealVar nsig("Nsig", "Nsig", 2349.0, 0, 2e5);

  RooRealVar width(Gamma(), Gamma(), 4.266, "MeV/c^{2}");
  RooRealVar mean(Mean(), Mean(), 1019.461, 1015.0, 1025.0, "MeV/c^{2}");
  RooRealVar sigma(Sigma(), Sigma(), 1.0, 0.05, 2.5, "MeV/c^{2}");
  RooRealVar q("q", "q", 1.5, "");

  auto& det = *bindPdf("det", qGaussian, x, RooConst(0), sigma, q);
  auto& bw = *bindPdf("bw", NA49RelativisticBW, x, mean, width);
  RooFFTConvPdf signal("signal", "signal", x, bw, det);
  RooExtendPdf esig("esig", "esig", signal, nsig, "norm");

  RooFormulaVar mX("mX", "@1 + @2 - @0", RooArgList(x, RooConst(xmax), RooConst(GetKThreshold())));
  RooRealVar k(GetPolyParName(0), GetPolyParName(0), -0.2, -10.0, 10.0);
  RooRealVar p(GetPolyParName(1), GetPolyParName(1), 0.2, 0.0, 1.0);
  RooArgusBG bkg("bkg", "bkg", mX, RooConst(xmax), k, p);

  RooRealVar nbkg("nbkg", "nbkg", 9384, 0, 2e5);
  RooExtendPdf ebkg("ebkg", "ebkg", bkg, nbkg);

  RooAddPdf model("model", "model", RooArgList(ebkg, esig));

  TStopwatch watch;

  RooFitResult* result = nullptr;
  for (int i = 0; i < 10; ++i)
    result = model.fitTo(data, PrintLevel(-1), Save(), Minimizer("Minuit2", "Minimize"));
  result->Print("V");

  // cout << int(watch.CpuTime() * 1000) << endl;
  cout << int(watch.RealTime() * 1000) << endl;
}