GammaFunctionGenerator.cc

Go to the documentation of this file.

//FAMOS headers
#include "FastSimulation/Utilities/interface/GammaFunctionGenerator.h"
#include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"

GammaFunctionGenerator::GammaFunctionGenerator() {
  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() {}

double GammaFunctionGenerator::shoot(RandomEngineAndDistribution const* random) const {
  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.;
  }
}

double GammaFunctionGenerator::gammaFrac(RandomEngineAndDistribution const* random) const {
  /* 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 {
  // 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) {
  //  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;
}