GammaContinuedFraction.cc File Reference

#include "CommonTools/Statistics/src/GammaContinuedFraction.h"
#include <cmath>
#include <iostream>

Go to the source code of this file.

Macros
#define	EPS   3.0e-7
#define	FPMIN   1.0e-30
#define	ITMAX   100

Functions
float	GammaContinuedFraction (float a, float x)

Macro Definition Documentation

#define EPS   3.0e-7

Definition at line 6 of file GammaContinuedFraction.cc.

Referenced by SingleParticleEvent::Eloss(), GammaContinuedFraction(), RectangularPixelTopology::localPosition(), reco::PreshowerCluster::operator==(), and SingleParticleEvent::subtractEloss().

#define FPMIN   1.0e-30

Definition at line 7 of file GammaContinuedFraction.cc.

Referenced by GammaContinuedFraction().

#define ITMAX   100

Definition at line 5 of file GammaContinuedFraction.cc.

Referenced by GammaContinuedFraction().

Function Documentation

float GammaContinuedFraction	(	float	a,
		float	x
	)

Returns the continued fraction summation of the (complement of the) incomplete gamma function P(a,x) evaluated by its continued fraction representation. source: Numerical Recipes

Definition at line 9 of file GammaContinuedFraction.cc.

References a, b, c, EcnaPython_AdcPeg12_S1_10_R170298_1_0_150_Dee0::cerr, ztail::d, EPS, FPMIN, h, mps_fire::i, and ITMAX.

Referenced by IncompleteGammaComplement::ln(), and IncompleteGammaComplement::value().

                                               {
  int i;
  float an, del;

  /* Set up for evaluating continued fraction by modified Lentz's method (par.5.2
     in Numerical Recipes in C) with b_0 = 0 */
  double b = x + 1.0 - a;
  double c = 1.0 / FPMIN;
  double d = 1.0 / b;
  double h = d;
  for (i = 1; i <= ITMAX; i++) {
    an = -i * (i - a);
    b += 2.0;
    d = an * d + b;
    if (fabs(d) < FPMIN)
      d = FPMIN;
    c = b + an / c;
    if (fabs(c) < FPMIN)
      c = FPMIN;
    d = 1.0 / d;
    del = d * c;
    h *= del;
    if (fabs(del - 1.0) < EPS)
      break;
  }
  if (i > ITMAX)
    std::cerr << "GammaContinuedFraction::a too large, "
              << "ITMAX too small" << std::endl;
  return h;
}