BCMath Namespace Reference

Some useful mathematic functions. More...

Functions
double	LogGaus (double x, double mean=0, double sigma=1, bool norm=false)
double	LogPoisson (double x, double par)
double	ApproxBinomial (int n, int k, double p)
double	LogApproxBinomial (int n, int k, double p)
double	LogBinomFactor (int n, int k)
double	ApproxLogFact (double x)
double	LogNoverK (int n, int k)
double	LogFact (int n)
int	Max (int x, int y)
double	Max (double x, double y)
int	Min (int x, int y)
double	Min (double x, double y)
int	Nint (double x)
double	rms (int n, const double *a)
double	LogBreitWignerNonRel (double x, double mean, double Gamma, bool norm=false)
double	LogBreitWignerRel (double x, double mean, double Gamma)
double	LogChi2 (double x, int n)
double	LogVoigtian (double x, double sigma, double gamma)
void	RandomChi2 (std::vector< double > &randoms, int K)
double	LogGaus (double x, double mean, double sigma, bool norm)
double	LogPoisson (double x, double par)
double	ApproxBinomial (int n, int k, double p)
double	LogApproxBinomial (int n, int k, double p)
double	LogBinomFactor (int n, int k)
double	ApproxLogFact (double x)
double	LogFact (int n)
double	LogNoverK (int n, int k)
int	Nint (double x)
double	rms (int n, const double *a)
double	LogBreitWignerNonRel (double x, double mean, double Gamma, bool norm)
double	LogBreitWignerRel (double x, double mean, double Gamma)
double	LogChi2 (double x, int n)
double	LogVoigtian (double x, double sigma, double gamma)
void	RandomChi2 (std::vector< double > &randoms, int K)

---

Detailed Description

Some useful mathematic functions.

Author:

Daniel Kollar

Kevin Kröninger

Jing Liu

Version:

1.0

Date:

08.2008 A namespace which encapsulates some mathematical functions necessary for BAT.

---

Function Documentation

double BCMath::ApproxBinomial	(	int	n,
		int	k,
		double	p
	)

Calculates Binomial probability using approximations for factorial calculations if calculation for number greater than 20 required using the BCMath::ApproxLogFact function.

Definition at line 62 of file BCMath.cxx.

{
   return exp( BCMath::LogApproxBinomial(n, k, p) );
}

double BCMath::ApproxBinomial	(	int	n,
		int	k,
		double	p
	)

Calculates Binomial probability using approximations for factorial calculations if calculation for number greater than 20 required using the BCMath::ApproxLogFact function.

Definition at line 62 of file BCMath.cxx.

{
   return exp( BCMath::LogApproxBinomial(n, k, p) );
}

double BCMath::ApproxLogFact

(

double

x

)

Calculates natural logarithm of the n-factorial (n!) using Srinivasa Ramanujan approximation log(n!) = n*log(n) - n + log(n*(1.+4.*n*(1.+2.*n)))/6. + log(PI)/2. if n > 20. If n <= 20 it uses BCMath::LogFact to calculate it exactly.

Definition at line 113 of file BCMath.cxx.

{
   if(x>BCMATH_NFACT_ALIMIT)
      return x*log(x) - x + log(x*(1.+4.*x*(1.+2.*x)))/6. + log(M_PI)/2.;

   else
      return BCMath::LogFact((int)x);
}

double BCMath::ApproxLogFact

(

double

x

)

Calculates natural logarithm of the n-factorial (n!) using Srinivasa Ramanujan approximation log(n!) = n*log(n) - n + log(n*(1.+4.*n*(1.+2.*n)))/6. + log(PI)/2. if n > 20. If n <= 20 it uses BCMath::LogFact to calculate it exactly.

Definition at line 113 of file BCMath.cxx.

{
   if(x>BCMATH_NFACT_ALIMIT)
      return x*log(x) - x + log(x*(1.+4.*x*(1.+2.*x)))/6. + log(M_PI)/2.;

   else
      return BCMath::LogFact((int)x);
}

double BCMath::LogApproxBinomial	(	int	n,
		int	k,
		double	p
	)

Calculates natural logarithm of the Binomial probability using approximations for factorial calculations if calculation for number greater than 20 required using the BCMath::ApproxLogFact function.

Definition at line 69 of file BCMath.cxx.

{
   // check p
   if (p == 0)
      return -1e99;

   else if (p == 1)
      return 0;

   // switch parameters if n < k
   if(n<k)
   {
      int a=n;
      n=k;
      k=a;
   }

   return BCMath::LogBinomFactor(n,k) + (double)k*log(p) + (double)(n-k)*log(1.-p);
}

double BCMath::LogApproxBinomial	(	int	n,
		int	k,
		double	p
	)

Calculates natural logarithm of the Binomial probability using approximations for factorial calculations if calculation for number greater than 20 required using the BCMath::ApproxLogFact function.

Definition at line 69 of file BCMath.cxx.

{
   // check p
   if (p == 0)
      return -1e99;

   else if (p == 1)
      return 0;

   // switch parameters if n < k
   if(n<k)
   {
      int a=n;
      n=k;
      k=a;
   }

   return BCMath::LogBinomFactor(n,k) + (double)k*log(p) + (double)(n-k)*log(1.-p);
}

double BCMath::LogBinomFactor	(	int	n,
		int	k
	)

Calculates natural logarithm of the Binomial factor "n over k" using approximations for factorial calculations if calculation for number greater than 20 required using the BCMath::ApproxLogFact function. Even for large numbers the calculation is performed precisely, if n-k < 5

Definition at line 91 of file BCMath.cxx.

{
   // switch parameters if n < k
   if(n<k)
   {
      int a=n;
      n=k;
      k=a;
   }

   if(n==k || k==0) return 0.;
   if(k==1 || k==n-1) return log((double)n);

   // if no approximation needed
   if(n<BCMATH_NFACT_ALIMIT || (n-k)<5) return BCMath::LogNoverK(n,k);

   // calculate final log(n over k) using approximations if necessary
   return BCMath::ApproxLogFact((double)n) - BCMath::ApproxLogFact((double)k) - BCMath::ApproxLogFact((double)(n-k));
}

double BCMath::LogBinomFactor	(	int	n,
		int	k
	)

Calculates natural logarithm of the Binomial factor "n over k" using approximations for factorial calculations if calculation for number greater than 20 required using the BCMath::ApproxLogFact function. Even for large numbers the calculation is performed precisely, if n-k < 5

Definition at line 91 of file BCMath.cxx.

{
   // switch parameters if n < k
   if(n<k)
   {
      int a=n;
      n=k;
      k=a;
   }

   if(n==k || k==0) return 0.;
   if(k==1 || k==n-1) return log((double)n);

   // if no approximation needed
   if(n<BCMATH_NFACT_ALIMIT || (n-k)<5) return BCMath::LogNoverK(n,k);

   // calculate final log(n over k) using approximations if necessary
   return BCMath::ApproxLogFact((double)n) - BCMath::ApproxLogFact((double)k) - BCMath::ApproxLogFact((double)(n-k));
}

double BCMath::LogBreitWignerNonRel	(	double	x,
		double	mean,
		double	Gamma,
		bool	norm = false
	)

Calculates the logarithm of the non-relativistic Breit-Wigner distribution.

Definition at line 209 of file BCMath.cxx.

{
   double bw = log(Gamma) - log( (x - mean) * (x - mean) + Gamma * Gamma / 4.0);

   if (norm)
      bw -= log(2. * M_PI);

   return bw;
}

double BCMath::LogBreitWignerNonRel	(	double	x,
		double	mean,
		double	Gamma,
		bool	norm = false
	)

Calculates the logarithm of the non-relativistic Breit-Wigner distribution.

Definition at line 209 of file BCMath.cxx.

{
   double bw = log(Gamma) - log( (x - mean) * (x - mean) + Gamma * Gamma / 4.0);

   if (norm)
      bw -= log(2. * M_PI);

   return bw;
}

double BCMath::LogBreitWignerRel	(	double	x,
		double	mean,
		double	Gamma
	)

Definition at line 221 of file BCMath.cxx.

{
   return -log( (x*x - mean*mean) * (x*x - mean*mean) + mean*mean*Gamma*Gamma );
}

double BCMath::LogBreitWignerRel	(	double	x,
		double	mean,
		double	Gamma
	)

Definition at line 221 of file BCMath.cxx.

{
   return -log( (x*x - mean*mean) * (x*x - mean*mean) + mean*mean*Gamma*Gamma );
}

double BCMath::LogChi2	(	double	x,
		int	n
	)

Calculates the logarithm of chi square function: chi2(double x; size_t n)

Definition at line 228 of file BCMath.cxx.

{
   if (x<0) {
      BCLog::OutWarning("BCMath::LogChi2 : parameter cannot be negative!");
      return -1e99;
   }

   if (x==0 && n==1) {
      BCLog::OutWarning("BCMath::LogChi2 : returned value is infinity!");
      return 1e99;
   }

   double nOver2 = ((double) n)/2.;

   return (nOver2-1.)*log(x) - x/2. - nOver2*log(2) - log(TMath::Gamma(nOver2));
}

double BCMath::LogChi2	(	double	x,
		int	n
	)

Calculates the logarithm of chi square function: chi2(double x; size_t n)

Definition at line 228 of file BCMath.cxx.

{
   if (x<0) {
      BCLog::OutWarning("BCMath::LogChi2 : parameter cannot be negative!");
      return -1e99;
   }

   if (x==0 && n==1) {
      BCLog::OutWarning("BCMath::LogChi2 : returned value is infinity!");
      return 1e99;
   }

   double nOver2 = ((double) n)/2.;

   return (nOver2-1.)*log(x) - x/2. - nOver2*log(2) - log(TMath::Gamma(nOver2));
}

double BCMath::LogFact

(

int

n

)

Calculates natural logarithm of the n-factorial (n!)

Definition at line 123 of file BCMath.cxx.

{
   double ln = 0.;

   for(int i=1;i<=n;i++)
      ln += log((double)i);

   return ln;
}

double BCMath::LogFact

(

int

n

)

Calculates natural logarithm of the n-factorial (n!)

Definition at line 123 of file BCMath.cxx.

{
   double ln = 0.;

   for(int i=1;i<=n;i++)
      ln += log((double)i);

   return ln;
}

double BCMath::LogGaus	(	double	x,
		double	mean = 0,
		double	sigma = 1,
		bool	norm = false
	)

Calculate the natural logarithm of a gaussian function with mean and sigma. If norm=true (default is false) the result is multiplied by the normalization constant, i.e. divided by sqrt(2*Pi)*sigma.

Definition at line 22 of file BCMath.cxx.

{
   // if we have a delta function, return fixed value
   if(sigma==0.)
      return 0;

   // if sigma is negative use absolute value
   if(sigma<0.)
      sigma *= -1.;

   double arg = (x-mean)/sigma;
   double result = -.5 * arg * arg;

   // check if we should add the normalization constant
   if(!norm)
      return result;

   // subtract the log of the denominator of the normalization constant
   // and return
   return result - log( sqrt( 2.* M_PI ) * sigma );
}

double BCMath::LogGaus	(	double	x,
		double	mean = 0,
		double	sigma = 1,
		bool	norm = false
	)

Calculate the natural logarithm of a gaussian function with mean and sigma. If norm=true (default is false) the result is multiplied by the normalization constant, i.e. divided by sqrt(2*Pi)*sigma.

Definition at line 22 of file BCMath.cxx.

{
   // if we have a delta function, return fixed value
   if(sigma==0.)
      return 0;

   // if sigma is negative use absolute value
   if(sigma<0.)
      sigma *= -1.;

   double arg = (x-mean)/sigma;
   double result = -.5 * arg * arg;

   // check if we should add the normalization constant
   if(!norm)
      return result;

   // subtract the log of the denominator of the normalization constant
   // and return
   return result - log( sqrt( 2.* M_PI ) * sigma );
}

double BCMath::LogNoverK	(	int	n,
		int	k
	)

Calculates natural logarithm of the Binomial factor "n over k".

Definition at line 135 of file BCMath.cxx.

{
   // switch parameters if n < k
   if(n<k)
   {
      int a = n;
      n = k;
      k = a;
   }

   if(n==k || k==0) return 0.;
   if(k==1 || k==n-1) return log((double)n);

   int lmax = BCMath::Max(k,n-k);
   int lmin = BCMath::Min(k,n-k);

   double ln = 0.;

   for(int i=n;i>lmax;i--)
      ln += log((double)i);
   ln -= BCMath::LogFact(lmin);

   return ln;
}

double BCMath::LogNoverK	(	int	n,
		int	k
	)

Calculates natural logarithm of the Binomial factor "n over k".

Definition at line 135 of file BCMath.cxx.

{
   // switch parameters if n < k
   if(n<k)
   {
      int a = n;
      n = k;
      k = a;
   }

   if(n==k || k==0) return 0.;
   if(k==1 || k==n-1) return log((double)n);

   int lmax = BCMath::Max(k,n-k);
   int lmin = BCMath::Min(k,n-k);

   double ln = 0.;

   for(int i=n;i>lmax;i--)
      ln += log((double)i);
   ln -= BCMath::LogFact(lmin);

   return ln;
}

double BCMath::LogPoisson	(	double	x,
		double	par
	)

Calculate the natural logarithm of a poisson distribution.

Definition at line 46 of file BCMath.cxx.

{
   if (par > 899)
      return BCMath::LogGaus(x,par,sqrt(par),true);

   if (x<0)
      return 0;

   if (x == 0.)
      return  -par;

   return x*log(par)-par-BCMath::ApproxLogFact(x);
}

double BCMath::LogPoisson	(	double	x,
		double	par
	)

Calculate the natural logarithm of a poisson distribution.

Definition at line 46 of file BCMath.cxx.

{
   if (par > 899)
      return BCMath::LogGaus(x,par,sqrt(par),true);

   if (x<0)
      return 0;

   if (x == 0.)
      return  -par;

   return x*log(par)-par-BCMath::ApproxLogFact(x);
}

double BCMath::LogVoigtian	(	double	x,
		double	sigma,
		double	gamma
	)

Calculates the logarithm of normalized voigtian function: voigtian(double x, double sigma, double gamma)

voigtian is a convolution of the following two functions: gaussian(x) = 1/(sqrt(2*pi)*sigma) * exp(x*x/(2*sigma*sigma) and lorentz(x) = (1/pi)*(gamma/2) / (x*x + (gamma/2)*(gamma/2))

it is singly peaked at x=0. The width of the peak is decided by sigma and gamma, so they should be positive.

Definition at line 247 of file BCMath.cxx.

{
   if (sigma<=0 || gamma<=0) {
      BCLog::OutWarning("BCMath::LogVoigtian : widths are negative or zero!");
      return -1e99;
   }

   return log(TMath::Voigt(x,sigma,gamma));
}

double BCMath::LogVoigtian	(	double	x,
		double	sigma,
		double	gamma
	)

Calculates the logarithm of normalized voigtian function: voigtian(double x, double sigma, double gamma)

voigtian is a convolution of the following two functions: gaussian(x) = 1/(sqrt(2*pi)*sigma) * exp(x*x/(2*sigma*sigma) and lorentz(x) = (1/pi)*(gamma/2) / (x*x + (gamma/2)*(gamma/2))

it is singly peaked at x=0. The width of the peak is decided by sigma and gamma, so they should be positive.

Definition at line 247 of file BCMath.cxx.

{
   if (sigma<=0 || gamma<=0) {
      BCLog::OutWarning("BCMath::LogVoigtian : widths are negative or zero!");
      return -1e99;
   }

   return log(TMath::Voigt(x,sigma,gamma));
}

double BCMath::Max	(	double	x,
		double	y
	)			[inline]

Definition at line 96 of file BCMath.h.

      { return x >= y ? x : y; }

int BCMath::Max	(	int	x,
		int	y
	)			[inline]

Returns the "greater or equal" of two numbers

Definition at line 93 of file BCMath.h.

      { return x >= y ? x : y; }

double BCMath::Min	(	double	x,
		double	y
	)			[inline]

Definition at line 105 of file BCMath.h.

   { return x <= y ? x : y; }

int BCMath::Min	(	int	x,
		int	y
	)			[inline]

Returns the "less or equal" of two numbers

Definition at line 102 of file BCMath.h.

      { return x <= y ? x : y; }

int BCMath::Nint

(

double

x

)

Returns the nearest integer of a double number.

Definition at line 162 of file BCMath.cxx.

{
   // round to integer
   int i;

   if (x >= 0)
   {
      i = (int)(x + .5);

      if (x + .5 == (double)i && i & 1)
         i--;
   }

   else
   {
      i = int(x - 0.5);

      if (x - 0.5 == double(i) && i & 1)
         i++;
   }

   return i;
}

int BCMath::Nint

(

double

x

)

Returns the nearest integer of a double number.

Definition at line 162 of file BCMath.cxx.

{
   // round to integer
   int i;

   if (x >= 0)
   {
      i = (int)(x + .5);

      if (x + .5 == (double)i && i & 1)
         i--;
   }

   else
   {
      i = int(x - 0.5);

      if (x - 0.5 == double(i) && i & 1)
         i++;
   }

   return i;
}

void BCMath::RandomChi2	(	std::vector< double > &	randoms,
		int	K
	)

Get N random numbers distributed according to chi square function with K degrees of freedom

Definition at line 264 of file BCMath.cxx.

                                                        {

   //fixed upper cutoff to 1000, might be too small
   TF1 *f = new TF1("chi2", chi2, 0.0, 1000, 1);
   f->SetParameter(0, K);
   f->SetNpx(500);
   //uses inverse-transform method
   //fortunately CDF only built once
   for (unsigned int i = 0; i < randoms.size(); i++)
      randoms.at(i) = f->GetRandom();
   delete f;
}

void BCMath::RandomChi2	(	std::vector< double > &	randoms,
		int	K
	)

Get N random numbers distributed according to chi square function with K degrees of freedom

Definition at line 264 of file BCMath.cxx.

                                                        {

   //fixed upper cutoff to 1000, might be too small
   TF1 *f = new TF1("chi2", chi2, 0.0, 1000, 1);
   f->SetParameter(0, K);
   f->SetNpx(500);
   //uses inverse-transform method
   //fortunately CDF only built once
   for (unsigned int i = 0; i < randoms.size(); i++)
      randoms.at(i) = f->GetRandom();
   delete f;
}

double BCMath::rms	(	int	n,
		const double *	a
	)

Returns the rms of an array.

Definition at line 188 of file BCMath.cxx.

{
   if (n <= 0 || !a)
      return 0;

   double sum = 0., sum2 = 0.;

   for (int i = 0; i < n; i++)
   {
      sum  += a[i];
      sum2 += a[i] * a[i];
   }

   double n1 = 1./(double)n;
   double mean = sum * n1;

   return sqrt(fabs(sum2 * n1 - mean * mean));
}

double BCMath::rms	(	int	n,
		const double *	a
	)

Returns the rms of an array.

Definition at line 188 of file BCMath.cxx.

{
   if (n <= 0 || !a)
      return 0;

   double sum = 0., sum2 = 0.;

   for (int i = 0; i < n; i++)
   {
      sum  += a[i];
      sum2 += a[i] * a[i];
   }

   double n1 = 1./(double)n;
   double mean = sum * n1;

   return sqrt(fabs(sum2 * n1 - mean * mean));
}

---