double negative_binomial_pdf(int k, double p, double r);

void test_NBD()
{
    
    std::random_device rd;
    std::mt19937 gen(rd());

	TH1F* hNBD = new TH1F("hNBD","hNBD",100,-0.5,100);
    TH1F* hNBDstd = new TH1F("hNBDstd","hNBDstd",100,-0.5,100);

	TF1* bin_rate_dist = new TF1("bin","ROOT::Math::negative_binomial_pdf(x[0],0.920519,18.9242)",0,100);
	bin_rate_dist->SetNpx(100);
    
    //std::negative_binomial_distribution
    
    std::negative_binomial_distribution<int> d(round(18.9242),0.920519);

	//mean: 1.63398

    //TF1* bin_rate_dist_v2 = new TF1("bin2","negative_binomial_pdf(x[0],0.920519,18.9242)",0,100);
	//bin_rate_dist_v2->SetNpx(1000);

	//Double_t par[2];
    //par[0] = 0.920519;
	//par[1] = 18.9242;
	//par[0] = 0.5;
	//par[1] = 3.0;

	//TF1* bin_rate_dist_v2 = new TF1("bin2","NBD",0.0,100.0,2);
	//bin_rate_dist_v2->SetNpx(100);
	//bin_rate_dist_v2->SetParameters(par);


	for(auto i=0;i<100000;++i){
		float rand = bin_rate_dist->GetRandom();
		//std::cout<<" Sample "<<rand<<"\n";
		hNBD->Fill(rand);
        hNBDstd->Fill(d(gen));
        
	}

	std::cout<<bin_rate_dist->Mean(0,100)<<"\n";
    std::cout<<" ROOT Histo mean "<<hNBD->GetMean()<<"\n";
    std::cout<<" STD Histo mean "<<hNBDstd->GetMean()<<"\n";
    

	TCanvas *c1 = new TCanvas();
	c1->Divide(2,2);
	c1->cd(1);
	hNBD->Draw();
	c1->cd(2);
	bin_rate_dist->Draw();
    c1->cd(3);
    hNBDstd->Draw();
    



}

double negative_binomial_pdf(int k, double p, double r) {

   if(r<0) return 0.0;
   if(p<0 || p>1.0) return 0.0;
   double a = ROOT::Math::lgamma(k+r) - ROOT::Math::lgamma(k + 1) - ROOT::Math::lgamma(r); 
   return std::exp(a) * std::pow(p, r) * std::pow(1-p, (double) k);
  
}

double NBD(Int_t *x,Double_t *par)
{
	if(par[1]<0) return 0.0;
        if(par[0]<0 || par[0]>1.0) return 0.0;
        double a = ROOT::Math::lgamma(x[0]+par[1]) - ROOT::Math::lgamma(x[0] + 1) - ROOT::Math::lgamma(par[1]); 
        return std::exp(a) * std::pow(par[0], par[1]) * std::pow(1-par[0], (double) x[0]);


}