#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, 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 31 of file GammaFunctionGenerator.cc.

31 {}

Member Function Documentation

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

private

values 0<a<1.

Definition at line 62 of file GammaFunctionGenerator.cc.

References create_public_lumi_plots::exp, RandomEngineAndDistribution::flatShoot(), frac, dqm-mbProfile::log, AlCaHLTBitMon_ParallelJobs::p, lumiQueryAPI::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 90 of file GammaFunctionGenerator.cc.

References approxLimit, coreCoeff, coreProba, RandomEngineAndDistribution::flatShoot(), dqm-mbProfile::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 113 of file GammaFunctionGenerator.cc.

References a, alpha, approxLimit, b, badRange, beta, coreProba, frac, integralToApproxLimit, myIncompleteGamma, na, theGammas, 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 33 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.;
     }
 }