#include <BaseNumericalRandomGenerator.h>

Inheritance diagram for BaseNumericalRandomGenerator:

Public Member Functions
	BaseNumericalRandomGenerator (double xmin=0., double xmax=1., int n=1000, int iter=6)

virtual double	function (double x)=0

double	generate (RandomEngineAndDistribution const *) const
	The random generation according to function() More...

double	generateExp (RandomEngineAndDistribution const *) const

double	generateLin (RandomEngineAndDistribution const *) const

void	initialize ()
	The initialization (numerical integarion, inversion) More...

bool	setSubInterval (double x1, double x2)
	To shoot in a given interval. More...

virtual	~BaseNumericalRandomGenerator ()
	Default destructor. More...

Protected Attributes
double	deltar

std::vector< double >	f

int	iter

int	n

double	rmin

std::vector< double >	sampling

double	xmax

double	xmin

Private Attributes
int	m

Detailed Description

Definition at line 27 of file BaseNumericalRandomGenerator.h.

Constructor & Destructor Documentation

◆ BaseNumericalRandomGenerator()

BaseNumericalRandomGenerator::BaseNumericalRandomGenerator	(	double	xmin = `0.`,
		double	xmax = `1.`,
		int	n = `1000`,
		int	iter = `6`
	)

Constructor that perform the necessary integration and inversion steps xmin and xmax are the generation bounds, n is the internal table size and iter is the number of iterations for the numerical part.

Definition at line 7 of file BaseNumericalRandomGenerator.cc.

     : xmin(xmin), xmax(xmax), n(n), iter(iter) {
   f.resize(n);
   sampling.resize(n);
 }

References f, n, and sampling.

◆ ~BaseNumericalRandomGenerator()

virtual BaseNumericalRandomGenerator::~BaseNumericalRandomGenerator ( )

inlinevirtual

Default destructor.

Definition at line 35 of file BaseNumericalRandomGenerator.h.

35 {}

Member Function Documentation

◆ function()

virtual double BaseNumericalRandomGenerator::function ( double x )

pure virtual

Implemented in HistogramGenerator, GammaNumericalGenerator, and LandauFluctuationGenerator.

◆ generate()

double BaseNumericalRandomGenerator::generate ( RandomEngineAndDistribution const * random ) const

The random generation according to function()

Definition at line 75 of file BaseNumericalRandomGenerator.cc.

                                                                                              {
   double r = rmin + deltar * random->flatShoot();
   int i = (int)r;
   double s = r - (double)i;
   //  cout << " i,r,s = " << i << " " << r << " " << s << endl;
   return sampling[i] + s * (sampling[i + 1] - sampling[i]);
 }

References RandomEngineAndDistribution::flatShoot(), mps_fire::i, createfilelist::int, alignCSCRings::r, rmin, alignCSCRings::s, and sampling.

Referenced by GammaNumericalGenerator::gamma(), and LandauFluctuationGenerator::landau().

◆ generateExp()

double BaseNumericalRandomGenerator::generateExp ( RandomEngineAndDistribution const * random ) const

The random generation according to function(), refined to generate as an exponential in each of the intervals

Definition at line 83 of file BaseNumericalRandomGenerator.cc.

                                                                                                 {
   double r = rmin + deltar * random->flatShoot();
   int i = (int)r;
   //  double s=r-(double)i;
   //  cout << " i,r,s = " << i << " " << r << " " << s << endl;
   double b = -std::log(f[i + 1] / f[i]) / (sampling[i + 1] - sampling[i]);
   double a = f[i] * std::exp(b * sampling[i]);
  
   double umin = -a / b * std::exp(-b * sampling[i]);
   double umax = -a / b * std::exp(-b * sampling[i + 1]);
   double u = (umax - umin) * random->flatShoot() + umin;
  
   return -std::log(-b / a * u) / b;
 }

References a, b, JetChargeProducer_cfi::exp, f, RandomEngineAndDistribution::flatShoot(), mps_fire::i, createfilelist::int, dqm-mbProfile::log, alignCSCRings::r, rmin, and sampling.

Referenced by GammaNumericalGenerator::gamma_exp().

◆ generateLin()

double BaseNumericalRandomGenerator::generateLin ( RandomEngineAndDistribution const * random ) const

The random generation according to function(), refined to generate as a linear function in each of the intervals

Definition at line 98 of file BaseNumericalRandomGenerator.cc.

                                                                                                 {
   double r = rmin + deltar * random->flatShoot();
   int i = (int)r;
   //  double s=r-(double)i;
   //  cout << " i,r,s = " << i << " " << r << " " << s << endl;
   double a = (f[i + 1] - f[i]) / (sampling[i + 1] - sampling[i]);
   double b = f[i] - a * sampling[i];
  
   double umin = a * sampling[i] * sampling[i] / 2. + b * sampling[i];
   double umax = a * sampling[i + 1] * sampling[i + 1] / 2. + b * sampling[i + 1];
   double u = (umax - umin) * random->flatShoot() + umin;
  
   return (-b + std::sqrt(b * b + 2. * a * u)) / a;
 }

References a, b, f, RandomEngineAndDistribution::flatShoot(), mps_fire::i, createfilelist::int, alignCSCRings::r, rmin, sampling, and mathSSE::sqrt().

Referenced by GammaNumericalGenerator::gamma_lin().

◆ initialize()

void BaseNumericalRandomGenerator::initialize ( )

The initialization (numerical integarion, inversion)

Definition at line 13 of file BaseNumericalRandomGenerator.cc.

                                               {
   m = n - 1;
   rmin = 0.;
   deltar = (double)m - rmin;
  
   std::vector<double> a, y, z, xnew;
   a.resize(n);
   y.resize(n);
   z.resize(n);
   xnew.resize(n);
  
   double sig1 = 0.;
  
   // Initialize sampling
   double du = (xmax - xmin) / (float)m;
   sampling[0] = xmin;
   for (int i = 1; i < n; ++i)
     sampling[i] = sampling[i - 1] + du;
  
   // Starting point for iterations
   for (int it = 0; it < iter; ++it) {
     // Calculate function values
     for (int i = 0; i < n; ++i)
       f[i] = function(sampling[i]);
  
     // Calculate bin areas
     for (int i = 0; i < m; ++i)
       a[i] = (sampling[i + 1] - sampling[i]) * (f[i + 1] + f[i]) / 2.;
  
     // Calculate cumulative spectrum Y values
     y[0] = 0.;
     for (int i = 1; i < n; ++i)
       y[i] = y[i - 1] + a[i - 1];
  
     // Put equidistant points on y scale
     double dz = y[n - 1] / (float)m;
     z[0] = 0;
     for (int i = 1; i < n; ++i)
       z[i] = z[i - 1] + dz;
  
     // Determine spacinf of Z points in between Y points
     // From this, determine new X values and finally replace old values
     xnew[0] = sampling[0];
     xnew[n - 1] = sampling[n - 1];
     int k = 0;
     for (int i = 1; i < m; ++i) {
       while (y[k + 1] < z[i])
         ++k;
       double r = (z[i] - y[k]) / (y[k + 1] - y[k]);
       xnew[i] = sampling[k] + (sampling[k + 1] - sampling[k]) * r;
     }
  
     for (int i = 0; i < n; ++i)
       sampling[i] = xnew[i];
  
     sig1 = sig1 + y[m];
     // std::cout << "BaseNumericalRandomGenerator::Iteration # " << it+1
     // << " Integral = " << sig1/(float)(it+1)
     // << std::endl;
   }
 }

References a, PVValHelper::dz, f, dqmMemoryStats::float, mps_fire::i, iter, dqmdumpme::k, m, n, alignCSCRings::r, rmin, sampling, xmax, xmin, y, and z.

Referenced by GammaNumericalGenerator::GammaNumericalGenerator(), HistogramGenerator::HistogramGenerator(), and LandauFluctuationGenerator::LandauFluctuationGenerator().

◆ setSubInterval()

bool BaseNumericalRandomGenerator::setSubInterval	(	double	x1,
		double	x2
	)

To shoot in a given interval.

Definition at line 113 of file BaseNumericalRandomGenerator.cc.

                                                                       {
   if (x1 < xmin || x2 > xmax)
     return false;
   if (x1 > x2) {
     double tmp = x1;
     x1 = x2;
     x2 = tmp;
   }
  
   unsigned ibin = 1;
   for (; ibin < (unsigned)n && x1 > sampling[ibin]; ++ibin)
     ;
   unsigned ibin1 = ibin;
   for (; ibin < (unsigned)n && x2 > sampling[ibin]; ++ibin)
     ;
   unsigned ibin2 = ibin;
  
   //  cout << sampling[ibin1-1] << " " << x1 << " " << sampling[ibin1] << endl;
   // cout << sampling[ibin2-1] << " " << x2 << " " << sampling[ibin2] << endl;
  
   rmin = ibin1 + (x1 - sampling[ibin1]) / (sampling[ibin1] - sampling[ibin1 - 1]);
   deltar = ibin2 + (x2 - sampling[ibin2]) / (sampling[ibin2] - sampling[ibin2 - 1]) - rmin;
   //  cout << rmin << " " << deltar << endl;
   return true;
 }