#include <GammaFunctionGenerator.h>

Public Member Functions
	GammaFunctionGenerator (const RandomEngine *engine)
	Constructor. More...

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

double	shoot () const

virtual	~GammaFunctionGenerator ()
	Destructor. More...

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

double	gammaInt () 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

const RandomEngine *	random

std::vector < GammaNumericalGenerator >	theGammas

double	xmax

double	xmin

Detailed Description

Definition at line 21 of file GammaFunctionGenerator.h.

Constructor & Destructor Documentation

GammaFunctionGenerator::GammaFunctionGenerator ( const RandomEngine * engine )

Constructor.

Definition at line 5 of file GammaFunctionGenerator.cc.

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

                                                                          :
   random(engine)
 {
 
   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(random,(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 33 of file GammaFunctionGenerator.cc.

33 {}

Member Function Documentation

double GammaFunctionGenerator::gammaFrac ( ) const

private

values 0<a<1.

Definition at line 64 of file GammaFunctionGenerator.cc.

References create_public_lumi_plots::exp, RandomEngine::flatShoot(), frac, create_public_lumi_plots::log, AlCaHLTBitMon_ParallelJobs::p, lumiQueryAPI::q, random, 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 ( ) const

private

integer values

Definition at line 92 of file GammaFunctionGenerator.cc.

References approxLimit, coreCoeff, coreProba, RandomEngine::flatShoot(), create_public_lumi_plots::log, na, random, 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();
     }
   //  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 115 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 ( ) const

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

Definition at line 35 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 ()/beta;
     }
       else if (na == 0)
     {
       return gammaFrac ()/beta;
     }
       else
     {
       double gi=gammaInt ();
       double gf=gammaFrac ();
       return (gi+gf)/beta;
     }
     }
   else
     {
       // an other generator has to be used in such a case
       return -1.;
     }
 }