Utilities/src/GammaFunctionGenerator.cc

0001 //FAMOS headers
0002 #include "FastSimulation/Utilities/interface/GammaFunctionGenerator.h"
0003 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0004
0005 GammaFunctionGenerator::GammaFunctionGenerator() {
0006   xmax = 30.;
0007
0008   for (unsigned i = 1; i <= 12; ++i) {
0009     // For all let's put the limit at 2.*(alpha-1)      (alpha-1 is the max of the dist)
0010     approxLimit.push_back(2 * ((double)i));
0011     myIncompleteGamma.a().setValue((double)i);
0012     integralToApproxLimit.push_back(myIncompleteGamma(approxLimit[i - 1]));
0013     theGammas.push_back(GammaNumericalGenerator((double)i, 1., 0, approxLimit[i - 1] + 1.));
0014   }
0015   coreCoeff.push_back(0.);  // alpha=1 not used
0016   coreCoeff.push_back(1. / 8.24659e-01);
0017   coreCoeff.push_back(1. / 7.55976e-01);
0018   coreCoeff.push_back(1. / 7.12570e-01);
0019   coreCoeff.push_back(1. / 6.79062e-01);
0020   coreCoeff.push_back(1. / 6.65496e-01);
0021   coreCoeff.push_back(1. / 6.48736e-01);
0022   coreCoeff.push_back(1. / 6.25185e-01);
0023   coreCoeff.push_back(1. / 6.09188e-01);
0024   coreCoeff.push_back(1. / 6.06221e-01);
0025   coreCoeff.push_back(1. / 6.05057e-01);
0026 }
0027
0028 GammaFunctionGenerator::~GammaFunctionGenerator() {}
0029
0030 double GammaFunctionGenerator::shoot(RandomEngineAndDistribution const* random) const {
0031   if (alpha < 0.)
0032     return -1.;
0033   if (badRange)
0034     return xmin / beta;
0035   if (alpha < 12) {
0036     if (alpha == na) {
0037       return gammaInt(random) / beta;
0038     } else if (na == 0) {
0039       return gammaFrac(random) / beta;
0040     } else {
0041       double gi = gammaInt(random);
0042       double gf = gammaFrac(random);
0043       return (gi + gf) / beta;
0044     }
0045   } else {
0046     // an other generator has to be used in such a case
0047     return -1.;
0048   }
0049 }
0050
0051 double GammaFunctionGenerator::gammaFrac(RandomEngineAndDistribution const* random) const {
0052   /* This is exercise 16 from Knuth; see page 135, and the solution is
0053      on page 551.  */
0054
0055   double p, q, x, u, v;
0056   p = M_E / (frac + M_E);
0057   do {
0058     u = random->flatShoot();
0059     v = random->flatShoot();
0060
0061     if (u < p) {
0062       x = exp((1 / frac) * log(v));
0063       q = exp(-x);
0064     } else {
0065       x = 1 - log(v);
0066       q = exp((frac - 1) * log(x));
0067     }
0068   } while (random->flatShoot() >= q);
0069
0070   return x;
0071 }
0072
0073 double GammaFunctionGenerator::gammaInt(RandomEngineAndDistribution const* random) const {
0074   // Exponential distribution : no approximation
0075   if (na == 1) {
0076     return xmin - log(random->flatShoot());
0077   }
0078
0079   unsigned gn = na - 1;
0080
0081   // are we sure to be in the tail
0082   if (coreProba == 0.)
0083     return xmin - coreCoeff[gn] * log(random->flatShoot());
0084
0085   // core-tail interval
0086   if (random->flatShoot() < coreProba) {
0087     return theGammas[gn].gamma_lin(random);
0088   }
0089   //  std::cout << " Tail called " << std::endl;
0090   return approxLimit[gn] - coreCoeff[gn] * log(random->flatShoot());
0091 }
0092
0093 void GammaFunctionGenerator::setParameters(double a, double b, double xm) {
0094   //  std::cout << "Setting parameters " << std::endl;
0095   alpha = a;
0096   beta = b;
0097   xmin = xm * beta;
0098   if (xm > xmax) {
0099     badRange = true;
0100     return;
0101   }
0102   badRange = false;
0103   na = 0;
0104
0105   if (alpha > 0. && alpha < 12)
0106     na = (unsigned)floor(alpha);
0107
0108   frac = alpha - na;
0109   // Now calculate the probability to shoot between approxLimit and xmax
0110   // The Incomplete gamma is normalized to 1
0111   if (na <= 1)
0112     return;
0113
0114   myIncompleteGamma.a().setValue((double)na);
0115
0116   unsigned gn = na - 1;
0117   //  std::cout << " na " << na << " xm " << xm << " beta " << beta << " xmin " << xmin << " approxLimit " << approxLimit[gn] << std::endl;
0118   if (xmin > approxLimit[gn]) {
0119     coreProba = 0.;
0120   } else {
0121     double tmp = (xmin != 0.) ? myIncompleteGamma(xmin) : 0.;
0122     coreProba = (integralToApproxLimit[gn] - tmp) / (1. - tmp);
0123     theGammas[gn].setSubInterval(xmin, approxLimit[gn]);
0124   }
0125   //  std::cout << " Proba " << coreProba << std::endl;
0126 }