#include <GammaFunctionGenerator.h>

Public Member Functions
	GammaFunctionGenerator ()
	Constructor. More...

void	setParameters (double a, double b, double xm)
	The parameters must be set before shooting. More...

double	shoot (RandomEngineAndDistribution const *) const

virtual	~GammaFunctionGenerator ()
	Destructor. More...

Private Member Functions
double	gammaFrac (RandomEngineAndDistribution const *) const
	values 0<a<1. More...

double	gammaInt (RandomEngineAndDistribution const *) const
	integer values More...

Private Attributes
double	alpha

std::vector< double >	approxLimit

bool	badRange

double	beta

std::vector< double >	coreCoeff

double	coreProba

double	frac

std::vector< double >	integralToApproxLimit

Genfun::IncompleteGamma	myIncompleteGamma

unsigned	na

std::vector < GammaNumericalGenerator >	theGammas

double	xmax

double	xmin

Detailed Description

Definition at line 21 of file GammaFunctionGenerator.h.

Constructor & Destructor Documentation

GammaFunctionGenerator::GammaFunctionGenerator ( )

Constructor.

Definition at line 5 of file GammaFunctionGenerator.cc.

References approxLimit, coreCoeff, alignCSCRings::e, mps_fire::i, integralToApproxLimit, myIncompleteGamma, theGammas, and xmax.

                                                {
   xmax = 30.;
 
   for (unsigned i = 1; i <= 12; ++i) {
     // For all let's put the limit at 2.*(alpha-1)      (alpha-1 is the max of the dist)
     approxLimit.push_back(2 * ((double)i));
     myIncompleteGamma.a().setValue((double)i);
     integralToApproxLimit.push_back(myIncompleteGamma(approxLimit[i - 1]));
     theGammas.push_back(GammaNumericalGenerator((double)i, 1., 0, approxLimit[i - 1] + 1.));
   }
   coreCoeff.push_back(0.);  // alpha=1 not used
   coreCoeff.push_back(1. / 8.24659e-01);
   coreCoeff.push_back(1. / 7.55976e-01);
   coreCoeff.push_back(1. / 7.12570e-01);
   coreCoeff.push_back(1. / 6.79062e-01);
   coreCoeff.push_back(1. / 6.65496e-01);
   coreCoeff.push_back(1. / 6.48736e-01);
   coreCoeff.push_back(1. / 6.25185e-01);
   coreCoeff.push_back(1. / 6.09188e-01);
   coreCoeff.push_back(1. / 6.06221e-01);
   coreCoeff.push_back(1. / 6.05057e-01);
 }

GammaFunctionGenerator::~GammaFunctionGenerator ( )

virtual

Destructor.

Definition at line 28 of file GammaFunctionGenerator.cc.

28 {}

Member Function Documentation

double GammaFunctionGenerator::gammaFrac ( RandomEngineAndDistribution const * random ) const

private

values 0<a<1.

Definition at line 51 of file GammaFunctionGenerator.cc.

References funct::exp(), RandomEngineAndDistribution::flatShoot(), frac, log, AlCaHLTBitMon_ParallelJobs::p, submitPVResolutionJobs::q, findQualityFiles::v, and x.

Referenced by shoot().

                                                                                         {
   /* This is exercise 16 from Knuth; see page 135, and the solution is
      on page 551.  */
 
   double p, q, x, u, v;
   p = M_E / (frac + M_E);
   do {
     u = random->flatShoot();
     v = random->flatShoot();
 
     if (u < p) {
       x = exp((1 / frac) * log(v));
       q = exp(-x);
     } else {
       x = 1 - log(v);
       q = exp((frac - 1) * log(x));
     }
   } while (random->flatShoot() >= q);
 
   return x;
 }

double GammaFunctionGenerator::gammaInt ( RandomEngineAndDistribution const * random ) const

private

integer values

Definition at line 73 of file GammaFunctionGenerator.cc.

References approxLimit, coreCoeff, coreProba, RandomEngineAndDistribution::flatShoot(), log, na, theGammas, and xmin.

Referenced by shoot().

                                                                                        {
   // Exponential distribution : no approximation
   if (na == 1) {
     return xmin - log(random->flatShoot());
   }
 
   unsigned gn = na - 1;
 
   // are we sure to be in the tail
   if (coreProba == 0.)
     return xmin - coreCoeff[gn] * log(random->flatShoot());
 
   // core-tail interval
   if (random->flatShoot() < coreProba) {
     return theGammas[gn].gamma_lin(random);
   }
   //  std::cout << " Tail called " << std::endl;
   return approxLimit[gn] - coreCoeff[gn] * log(random->flatShoot());
 }

void GammaFunctionGenerator::setParameters	(	double	a,
		double	b,
		double	xm
	)

The parameters must be set before shooting.

Definition at line 93 of file GammaFunctionGenerator.cc.

References a, alpha, approxLimit, b, badRange, beta, coreProba, frac, integralToApproxLimit, myIncompleteGamma, na, theGammas, createJobs::tmp, xmax, and xmin.

Referenced by EMShower::compute().

                                                                         {
   //  std::cout << "Setting parameters " << std::endl;
   alpha = a;
   beta = b;
   xmin = xm * beta;
   if (xm > xmax) {
     badRange = true;
     return;
   }
   badRange = false;
   na = 0;
 
   if (alpha > 0. && alpha < 12)
     na = (unsigned)floor(alpha);
 
   frac = alpha - na;
   // Now calculate the probability to shoot between approxLimit and xmax
   // The Incomplete gamma is normalized to 1
   if (na <= 1)
     return;
 
   myIncompleteGamma.a().setValue((double)na);
 
   unsigned gn = na - 1;
   //  std::cout << " na " << na << " xm " << xm << " beta " << beta << " xmin " << xmin << " approxLimit " << approxLimit[gn] << std::endl;
   if (xmin > approxLimit[gn]) {
     coreProba = 0.;
   } else {
     double tmp = (xmin != 0.) ? myIncompleteGamma(xmin) : 0.;
     coreProba = (integralToApproxLimit[gn] - tmp) / (1. - tmp);
     theGammas[gn].setSubInterval(xmin, approxLimit[gn]);
   }
   //  std::cout << " Proba " << coreProba << std::endl;
 }

double GammaFunctionGenerator::shoot ( RandomEngineAndDistribution const * random ) const

shoot along a gamma distribution with shape parameter alpha and scale beta values > xmin

Definition at line 30 of file GammaFunctionGenerator.cc.

References alpha, badRange, beta, gammaFrac(), gammaInt(), gf, na, and xmin.

Referenced by EMShower::compute().

                                                                                     {
   if (alpha < 0.)
     return -1.;
   if (badRange)
     return xmin / beta;
   if (alpha < 12) {
     if (alpha == na) {
       return gammaInt(random) / beta;
     } else if (na == 0) {
       return gammaFrac(random) / beta;
     } else {
       double gi = gammaInt(random);
       double gf = gammaFrac(random);
       return (gi + gf) / beta;
     }
   } else {
     // an other generator has to be used in such a case
     return -1.;
   }
 }