Utilities/src/BaseNumericalRandomGenerator.cc

0001 #include "FastSimulation/Utilities/interface/BaseNumericalRandomGenerator.h"
0002 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0003
0004 #include <cmath>
0005 // #include <iostream>
0006
0007 BaseNumericalRandomGenerator::BaseNumericalRandomGenerator(double xmin, double xmax, int n, int iter)
0008     : xmin(xmin), xmax(xmax), n(n), iter(iter) {
0009   f.resize(n);
0010   sampling.resize(n);
0011 }
0012
0013 void BaseNumericalRandomGenerator::initialize() {
0014   m = n - 1;
0015   rmin = 0.;
0016   deltar = (double)m - rmin;
0017
0018   std::vector<double> a, y, z, xnew;
0019   a.resize(n);
0020   y.resize(n);
0021   z.resize(n);
0022   xnew.resize(n);
0023
0024   double sig1 = 0.;
0025
0026   // Initialize sampling
0027   double du = (xmax - xmin) / (float)m;
0028   sampling[0] = xmin;
0029   for (int i = 1; i < n; ++i)
0030     sampling[i] = sampling[i - 1] + du;
0031
0032   // Starting point for iterations
0033   for (int it = 0; it < iter; ++it) {
0034     // Calculate function values
0035     for (int i = 0; i < n; ++i)
0036       f[i] = function(sampling[i]);
0037
0038     // Calculate bin areas
0039     for (int i = 0; i < m; ++i)
0040       a[i] = (sampling[i + 1] - sampling[i]) * (f[i + 1] + f[i]) / 2.;
0041
0042     // Calculate cumulative spectrum Y values
0043     y[0] = 0.;
0044     for (int i = 1; i < n; ++i)
0045       y[i] = y[i - 1] + a[i - 1];
0046
0047     // Put equidistant points on y scale
0048     double dz = y[n - 1] / (float)m;
0049     z[0] = 0;
0050     for (int i = 1; i < n; ++i)
0051       z[i] = z[i - 1] + dz;
0052
0053     // Determine spacinf of Z points in between Y points
0054     // From this, determine new X values and finally replace old values
0055     xnew[0] = sampling[0];
0056     xnew[n - 1] = sampling[n - 1];
0057     int k = 0;
0058     for (int i = 1; i < m; ++i) {
0059       while (y[k + 1] < z[i])
0060         ++k;
0061       double r = (z[i] - y[k]) / (y[k + 1] - y[k]);
0062       xnew[i] = sampling[k] + (sampling[k + 1] - sampling[k]) * r;
0063     }
0064
0065     for (int i = 0; i < n; ++i)
0066       sampling[i] = xnew[i];
0067
0068     sig1 = sig1 + y[m];
0069     // std::cout << "BaseNumericalRandomGenerator::Iteration # " << it+1
0070     // << " Integral = " << sig1/(float)(it+1)
0071     // << std::endl;
0072   }
0073 }
0074
0075 double BaseNumericalRandomGenerator::generate(RandomEngineAndDistribution const* random) const {
0076   double r = rmin + deltar * random->flatShoot();
0077   int i = (int)r;
0078   double s = r - (double)i;
0079   //  cout << " i,r,s = " << i << " " << r << " " << s << endl;
0080   return sampling[i] + s * (sampling[i + 1] - sampling[i]);
0081 }
0082
0083 double BaseNumericalRandomGenerator::generateExp(RandomEngineAndDistribution const* random) const {
0084   double r = rmin + deltar * random->flatShoot();
0085   int i = (int)r;
0086   //  double s=r-(double)i;
0087   //  cout << " i,r,s = " << i << " " << r << " " << s << endl;
0088   double b = -std::log(f[i + 1] / f[i]) / (sampling[i + 1] - sampling[i]);
0089   double a = f[i] * std::exp(b * sampling[i]);
0090
0091   double umin = -a / b * std::exp(-b * sampling[i]);
0092   double umax = -a / b * std::exp(-b * sampling[i + 1]);
0093   double u = (umax - umin) * random->flatShoot() + umin;
0094
0095   return -std::log(-b / a * u) / b;
0096 }
0097
0098 double BaseNumericalRandomGenerator::generateLin(RandomEngineAndDistribution const* random) const {
0099   double r = rmin + deltar * random->flatShoot();
0100   int i = (int)r;
0101   //  double s=r-(double)i;
0102   //  cout << " i,r,s = " << i << " " << r << " " << s << endl;
0103   double a = (f[i + 1] - f[i]) / (sampling[i + 1] - sampling[i]);
0104   double b = f[i] - a * sampling[i];
0105
0106   double umin = a * sampling[i] * sampling[i] / 2. + b * sampling[i];
0107   double umax = a * sampling[i + 1] * sampling[i + 1] / 2. + b * sampling[i + 1];
0108   double u = (umax - umin) * random->flatShoot() + umin;
0109
0110   return (-b + std::sqrt(b * b + 2. * a * u)) / a;
0111 }
0112
0113 bool BaseNumericalRandomGenerator::setSubInterval(double x1, double x2) {
0114   if (x1 < xmin || x2 > xmax)
0115     return false;
0116   if (x1 > x2) {
0117     double tmp = x1;
0118     x1 = x2;
0119     x2 = tmp;
0120   }
0121
0122   unsigned ibin = 1;
0123   for (; ibin < (unsigned)n && x1 > sampling[ibin]; ++ibin)
0124     ;
0125   unsigned ibin1 = ibin;
0126   for (; ibin < (unsigned)n && x2 > sampling[ibin]; ++ibin)
0127     ;
0128   unsigned ibin2 = ibin;
0129
0130   //  cout << sampling[ibin1-1] << " " << x1 << " " << sampling[ibin1] << endl;
0131   // cout << sampling[ibin2-1] << " " << x2 << " " << sampling[ibin2] << endl;
0132
0133   rmin = ibin1 + (x1 - sampling[ibin1]) / (sampling[ibin1] - sampling[ibin1 - 1]);
0134   deltar = ibin2 + (x2 - sampling[ibin2]) / (sampling[ibin2] - sampling[ibin2 - 1]) - rmin;
0135   //  cout << rmin << " " << deltar << endl;
0136   return true;
0137 }