#include <GammaFunctionGenerator.h>

Public Member Functions
	GammaFunctionGenerator (const RandomEngine *engine)
	Constructor.
void	setParameters (double a, double b, double xm)
	The parameters must be set before shooting.
double	shoot () const
virtual	~GammaFunctionGenerator ()
	Destructor.
Private Member Functions
double	gammaFrac () const
	values 0<a<1.
double	gammaInt () const
	integer values
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, ExpressReco_HICollisions_FallBack::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.

{}

Member Function Documentation

double GammaFunctionGenerator::gammaFrac ( ) const [private]

values 0<a<1.

Definition at line 64 of file GammaFunctionGenerator.cc.

References funct::exp(), RandomEngine::flatShoot(), frac, funct::log(), L1TEmulatorMonitor_cff::p, lumiQueryAPI::q, random, v, and ExpressReco_HICollisions_FallBack::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(), funct::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.;
    }
}