File indexing completed on 2023-03-17 11:01:02
0001
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
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.);
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
0047 return -1.;
0048 }
0049 }
0050
0051 double GammaFunctionGenerator::gammaFrac(RandomEngineAndDistribution const* random) const {
0052
0053
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
0075 if (na == 1) {
0076 return xmin - log(random->flatShoot());
0077 }
0078
0079 unsigned gn = na - 1;
0080
0081
0082 if (coreProba == 0.)
0083 return xmin - coreCoeff[gn] * log(random->flatShoot());
0084
0085
0086 if (random->flatShoot() < coreProba) {
0087 return theGammas[gn].gamma_lin(random);
0088 }
0089
0090 return approxLimit[gn] - coreCoeff[gn] * log(random->flatShoot());
0091 }
0092
0093 void GammaFunctionGenerator::setParameters(double a, double b, double xm) {
0094
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
0110
0111 if (na <= 1)
0112 return;
0113
0114 myIncompleteGamma.a().setValue((double)na);
0115
0116 unsigned gn = na - 1;
0117
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
0126 }