#include <BaseNumericalRandomGenerator.h>

Inheritance diagram for BaseNumericalRandomGenerator:

Public Member Functions
	BaseNumericalRandomGenerator (const RandomEngine *engine, double xmin=0., double xmax=1., int n=1000, int iter=6)
virtual double	function (double x)=0
double	generate () const
	The random generation according to function()
double	generateExp () const
double	generateLin () const
void	initialize ()
	The initialization (numerical integarion, inversion)
bool	setSubInterval (double x1, double x2)
	To shoot in a given interval.
virtual	~BaseNumericalRandomGenerator ()
	Default destructor.
Protected Attributes
double	deltar
std::vector< double >	f
int	iter
int	n
const RandomEngine *	random
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	(	const RandomEngine *	engine,
		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.

References f, and sampling.

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

virtual BaseNumericalRandomGenerator::~BaseNumericalRandomGenerator ( ) [inline, virtual]

Default destructor.

Definition at line 41 of file BaseNumericalRandomGenerator.h.

{}

Member Function Documentation

virtual double BaseNumericalRandomGenerator::function ( double x ) [pure virtual]

Implemented in GammaNumericalGenerator, HistogramGenerator, and LandauFluctuationGenerator.

Referenced by initialize().

double BaseNumericalRandomGenerator::generate ( ) const

The random generation according to function()

Definition at line 84 of file BaseNumericalRandomGenerator.cc.

References deltar, RandomEngine::flatShoot(), i, alignCSCRings::r, random, rmin, alignCSCRings::s, and sampling.

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

                                             {

  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]);

}

double BaseNumericalRandomGenerator::generateExp ( ) const

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

Definition at line 95 of file BaseNumericalRandomGenerator.cc.

References a, b, deltar, funct::exp(), f, RandomEngine::flatShoot(), i, funct::log(), alignCSCRings::r, random, rmin, and sampling.

Referenced by GammaNumericalGenerator::gamma_exp().

                                                {

  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;

}

double BaseNumericalRandomGenerator::generateLin ( ) const

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

Definition at line 113 of file BaseNumericalRandomGenerator.cc.

References a, b, deltar, f, RandomEngine::flatShoot(), i, alignCSCRings::r, random, rmin, sampling, and mathSSE::sqrt().

Referenced by GammaNumericalGenerator::gamma_lin().

                                                {

  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;

}

void BaseNumericalRandomGenerator::initialize ( )

The initialization (numerical integarion, inversion)

Definition at line 18 of file BaseNumericalRandomGenerator.cc.

References a, deltar, f, function(), i, iter, gen::k, m, n, alignCSCRings::r, rmin, sampling, xmax, xmin, detailsBasic3DVector::y, and z.

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

                                         {
  
  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;

  }

}

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

To shoot in a given interval.

Definition at line 131 of file BaseNumericalRandomGenerator.cc.

References deltar, rmin, sampling, tmp, and xmax.

                                                                    {
  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;
}