#include "TCanvas.h"
#include "TH1D.h"
#include "TF1.h"
#include "Math/PdfFuncMathCore.h"

/*
This is the function:
ref https://root.cern/doc/master/group__PdfFunc.html#ga4a9cab35372227aa21c2de11f625325e

double negative_binomial_pdf(unsigned int k, double p, double n) {
	// impelment in term of gamma function
	 
	if (n < 0)  return 0.0;
	if (p < 0 || p > 1.0) return 0.0;
   
	double coeff = ROOT::Math::lgamma(k+n) - ROOT::Math::lgamma(k+1.0) - ROOT::Math::lgamma(n);
	return std::exp(coeff + n * std::log(p) + double(k) * ROOT::Math::log1p(-p));

}
*/

double my_negative_binomial_pdf(double* x, double* par)
{
	const unsigned int k = static_cast<unsigned int>(x[0]);
	const double p = par[0];
	const double n = par[1];
//	std::cout << "k " << k << ", p " << p << ", n " << n << std::endl;
	Double_t f = ROOT::Math::negative_binomial_pdf(k, p, n);
	return f;
}

void neg_binomial_plot(unsigned int runs=1000)
{
	TCanvas *c = new TCanvas("c", "negative binomial", 900, 700);
	TH1D* h = new TH1D("h", "negative binomial", 150, 0., 75.);
	
	const double XMin = 0.;
	const double XMax = 100.;
	const int NParm = 2;
	
	// cf. https://root.cern.ch/doc/master/classTF1.html#F4
	TF1* fNegBinomial = new TF1("fNegBinomial", my_negative_binomial_pdf, XMin, XMax, NParm);
	fNegBinomial->SetParameter(0, 0.25); // p
	fNegBinomial->SetParameter(1, 2000); // n
	
	for (unsigned int r=0; r<runs; r++)
	{	
		const double val = fNegBinomial->GetRandom();
		h->Fill(val);
	}
	
	c->cd();
	h->Draw();

}