#define N_BINS 4800
#define pshift 0.22278298

#ifndef FUNCTION_H
#define FUNCTION_H





float binomial_cal ( int nshot, int nexp )
{
	float p = 1.0/25.0;
	float factpart = 1;
	for ( int i=0; i<nexp; i++ )
		factpart = factpart * float( nshot-i ) / float( i+1 );
	float otherpart = pow( p, nexp ) * pow ( 1-p, nshot-nexp );
	float result = factpart * otherpart;
	return result;
}


double poisson_cal ( int nexp, double lambda )
{
	double result = pow( lambda, nexp ) * exp( -lambda );
	for ( int i=1; i<=nexp; i++ )
		result = result / i;
	return result;
}




double fval_current[N_BINS];
double fval_previous[N_BINS];
double fval_landau[N_BINS];
double fval_gauss[N_BINS];

double value;
double tau;
double fforward;
double fbackward;
double max_var;
int    n_lsteps;
int    n_gsteps;
double stepsize;

int  bin_zero;
int  bin_max;
int  bin_return;
int  n_stepstobound;
int  bin_topbound;


double fastconvolution (
		double var,
		int n_conv,
		double lmean,
		double lwidth,
		double gsigma )
{
	max_var = max ( lmean*7.0, 150.0 );
	bin_zero = 200;

	n_lsteps = min ( int( 4*max_var ), N_BINS-bin_zero );
	n_gsteps = min ( int( 20*gsigma ), N_BINS-bin_zero );
	stepsize = max_var/float(n_lsteps);

	bin_max  = bin_zero + n_lsteps;

	std::fill ( std::begin( fval_landau ), std::end( fval_landau ), 0.0 );
	std::fill ( std::begin( fval_gauss ), std::end( fval_gauss ), 0.0 );
	std::fill ( std::begin( fval_current ), std::end( fval_current ), 0.0 );
	std::fill ( std::begin( fval_previous ), std::end( fval_previous ), 0.0 );


	for ( int ibin=0; ibin<bin_max; ibin++ )
	{
		fval_gauss[ibin] = TMath::Gaus(
				(ibin+0.5-bin_zero)*stepsize, 0, gsigma, true );
		fval_landau[ibin] =  TMath::Landau(
				(ibin+0.5-bin_zero)*stepsize, lmean, lwidth, true );
		fval_previous[ibin] = TMath::Landau(
				(ibin+0.5-bin_zero)*stepsize, lmean, lwidth, true );
		fval_current[ibin] = TMath::Landau(
				(ibin+0.5-bin_zero)*stepsize, lmean, lwidth, true );
	}


	for ( int iconv=1; iconv<n_conv; iconv++ )
	{
		for ( int ibin=bin_zero; ibin!=bin_max; ibin++ )
		{
			n_stepstobound = ibin - bin_zero;
			bin_topbound = bin_zero + n_stepstobound/2;

			value = 0;
			for ( int step=bin_zero; step<=bin_topbound; step++ )
				value += (
						fval_previous[step-1]*fval_landau[ibin-step+bin_zero] +
						fval_previous[ibin-step+bin_zero]*fval_landau[step-1] );
			value += ((n_stepstobound)%2) *
				fval_previous[bin_zero+(n_stepstobound+1)/2] *
				fval_landau[bin_zero+(n_stepstobound+1)/2];

			fval_current[ibin] = value*stepsize;
		}

		memmove ( fval_previous, fval_current, N_BINS*sizeof(double) );
	}


	bin_return = int(var/stepsize)+bin_zero;
	for ( int ibin=bin_zero; ibin<bin_max; ibin++ )
	{
		value = !n_gsteps*fval_previous[ibin]/stepsize;
		for ( int step=-n_gsteps; step<n_gsteps; step++ )
			value += ( fval_gauss[step+bin_zero+1]*fval_previous[ibin-step-1] );

		fval_current[ibin] = value*stepsize;
	}

	return 1.0*fval_current[bin_return];
}


double simpleconvolution ( double var, double lmean, double lwidth )
{
	double bound1 = 0.0;
	max_var = lmean*6;
	stepsize = max_var/200;
	double result = 0;

	for ( int i=0; i<200; i++ )
	{
		tau = ( i+0.5 ) * stepsize;
		fforward  = TMath::Landau( tau, lmean, lwidth, true );
		fbackward = TMath::Landau( var-tau, lmean, lwidth, true );
		value += fforward*fbackward;
	}

	int varbin = var / stepsize;

	return fval_current[varbin];
}


double pdfsum_multilandau ( double *var, double *par )
{
	double adc = var[0];
	double lmean  = par[0]+pshift*par[1];
	double lwidth = par[1];
	double area   = par[2];
	double flux   = par[3];
	double sigped = par[4];
	double sigphy = par[5];
	double result;

	//if ( lwidth<0 || lmean<lwidth )   return 0;

	double sigma1 = sqrt ( sigped*sigped +   sigphy*sigphy );
	double sigma2 = sqrt ( sigped*sigped + 2*sigphy*sigphy );
	double sigma3 = sqrt ( sigped*sigped + 3*sigphy*sigphy );
	//double sigma4 = sqrt ( sigped*sigped + 4*sigphy*sigphy );

	double flandau1 = fastconvolution ( adc, 1, lmean, lwidth, sigma1 );
	double flandau2 = fastconvolution ( adc, 2, lmean, lwidth, sigma2 );
	double flandau3 = fastconvolution ( adc, 3, lmean, lwidth, sigma3 );
	//double flandau4 = fastconvolution ( adc, 4, lmean, lwidth, sigma4 );

	//double chance1 = binomial_cal ( flux, 1 );
	//double chance2 = binomial_cal ( flux, 2 );
	//double chance3 = binomial_cal ( flux, 3 );
	//double chance4 = binomial_cal ( flux, 4 );
	double chance1 = ROOT::Math::binomial_pdf( 1, 0.04, flux );
	double chance2 = ROOT::Math::binomial_pdf( 2, 0.04, flux );
	double chance3 = ROOT::Math::binomial_pdf( 3, 0.04, flux );

	//result = area * (
			//chance1*flandau1 +
			//chance2*flandau2 +
			//chance3*flandau3 +
			//chance4*flandau4 );
	result = area * (
			chance1*flandau1 +
			chance2*flandau2 +
			chance3*flandau3 );
	return result;
}


double pdfsum_doublelandau ( double *var, double *par )
{
	int n_landau_step = 100;
	int n_gauss_step  = 100;
	
	double adc = var[0];
	double lmean  = par[0]+pshift*par[1];
	double lwidth = par[1];
	double area   = par[2];
	double flux   = par[3];
	double sigped = par[4];
	double sigphy = par[5];
	double result;

	double sigma1 = sqrt( sigped*sigped +   sigphy*sigphy );
	double sigma2 = sqrt( sigped*sigped + 2*sigphy*sigphy );
	//double chance1 = binomial_cal( flux, 1 );
	//double chance2 = binomial_cal( flux, 2 );
	double chance1 = ROOT::Math::binomial_pdf( 1, 0.04, flux );
	double chance2 = ROOT::Math::binomial_pdf( 2, 0.04, flux );

	double fconv1hit = 0.0;
	double lower_adc1 = adc - 5*sigma1;
	double g_stepsize1 = 10*sigma1 / n_gauss_step;
	for ( int gstep=0; gstep<n_gauss_step; gstep++ )
	{
		double ladc = lower_adc1 + ( gstep+0.5 )*g_stepsize1;
		double gadc = adc - ladc;

		double flandau = TMath::Landau( ladc, lmean, lwidth, true );
		double fgauss  = TMath::Gaus( gadc, 0, sigma1, true );
		fconv1hit += ( flandau*fgauss );
	}

	double fconv2hit = 0.0;
	double lower_adc2 = adc - 5*sigma2;
	double g_stepsize2 = 10*sigma2 / n_gauss_step;
	for ( int gstep=0; gstep<n_gauss_step; gstep++ )
	{
		double ladc = lower_adc2 + ( gstep+0.5 )*g_stepsize2;
		double gadc = adc - ladc;

		double fconv2landau = 0.0;
		double l_stepsize = 8*lmean / n_landau_step;
		int    max_landau_step = ladc / l_stepsize;
		for ( int lstep=0; lstep<max_landau_step; lstep++ )
		{
			double subadc1 = ( lstep+0.5 )*l_stepsize;
			double subadc2 = ladc - subadc1;

			double flandau1 = TMath::Landau( subadc1, lmean, lwidth, true );
			double flandau2 = TMath::Landau( subadc2, lmean, lwidth, true );
			fconv2landau += ( flandau1*flandau2 );
		}
		double fgauss = TMath::Gaus( gadc, 0, sigma2, true );

		fconv2hit += ( fconv2landau*l_stepsize*fgauss );
	}

	result = area * (
			chance1*fconv1hit*g_stepsize1 +
			chance2*fconv2hit*g_stepsize2 );
	return result;
}


#endif