Heed::PointsRan Class Reference

#include <PointsRan.h>

Public Member Functions
	PointsRan ()
	PointsRan (const std::vector< double > &fx, const std::vector< double > &fy, double fxmin, double fxmax)
double	ran (double flat_ran) const
double	get_integ_total () const
double	get_integ_active () const
void	print (std::ostream &file) const

Detailed Description

Generates random numbers according to array. This class generates random numbers according to a pointwise distribution, with linear interpolation between given points and also with linear extrapolation.

Definition at line 24 of file PointsRan.h.

Constructor & Destructor Documentation

PointsRan() [1/2]

Heed::PointsRan::PointsRan ( )

inline

Definition at line 46 of file PointsRan.h.

{}

PointsRan() [2/2]

Heed::PointsRan::PointsRan	(	const std::vector< double > &	fx,
		const std::vector< double > &	fy,
		double	fxmin,
		double	fxmax
	)

Constructor

Parameters

fx	x values
fy	y values
fxmin	Minimum of generated distribution. If less then x[0], extend the distribution by linear extrapolation. If the extrapolated line crosses zero, the extension stops there. Otherwise it stops at fxmin.
fxmax	Maximum of generated distribution. If greater than x[q-1], extend the distribution by linear extrapolation.

Definition at line 19 of file PointsRan.cpp.

    : xmin(fxmin), xmax(fxmax), x(fx), y(fy) {
  mfunnamep("PointsRan::PointsRan(...)");
  check_econd12(x.size(), !=, y.size(), mcerr);
  check_econd11(x.size(), < 2, mcerr);
  check_econd12(xmin, >=, xmax, mcerr);
  const long q = x.size();
  for (long n = 0; n < q - 1; n++) {
    check_econd12(x[n], >=, x[n + 1], mcerr);
  }
  for (long n = 0; n < q; n++) {
    check_econd11(y[n], < 0.0, mcerr);
  }
  iy.resize(q);
  a.resize(q - 1);
  for (long n = 0; n < q - 1; n++) {
    a[n] = (y[n + 1] - y[n]) / (x[n + 1] - x[n]);
  }
  if (xmin < x[0]) {
    double y0 = y[0] - a[0] * (x[0] - xmin);
    if (y0 < 0.0) {
      x[0] = x[0] - y[0] / a[0];
      check_econd12(xmax, <, x[0], mcerr);
      xmin = x[0];
      y[0] = 0.0;
    } else {
      y[0] = y0;
      x[0] = xmin;
    }
  }
  if (xmax > x[q - 1]) {
    double yq = y[q - 1] + a[q - 2] * (xmax - x[q - 1]);
    if (yq < 0.0) {
      x[q - 1] = x[q - 1] - y[q - 1] / a[0];
      check_econd12(xmin, >, x[q - 1], mcerr);
      xmax = x[q - 1];
      y[q - 1] = 0.0;
    } else {
      x[q - 1] = xmax;
      y[q - 1] = yq;
    }
  }
  // now xmin and xmax should be between x[0] and x[q-1]
  iy[0] = 0.0;
  integ_start = 0.0;
  n_start = 0;
  for (long n = 1; n < q; n++) {
    iy[n] = iy[n - 1] + 0.5 * (x[n] - x[n - 1]) * (y[n - 1] + y[n]);
    if (xmin >= x[n]) {
      integ_start = iy[n];
      n_start = n;
    } else if (xmin > x[n - 1]) {
      integ_start +=
          (xmin - x[n - 1]) * (y[n - 1] + 0.5 * a[n - 1] * (xmin - x[n - 1]));
      n_start = n - 1;
    }
    if (xmax >= x[n]) {
      integ_finish = iy[n];
      n_finish = n;
    } else if (xmax > x[n - 1]) {
      integ_finish +=
          (xmax - x[n - 1]) * (y[n - 1] + 0.5 * a[n - 1] * (xmax - x[n - 1]));
      n_finish = n;
    }
  }
  integ_active = integ_finish - integ_start;
  double s = iy[q - 1];
  integ_total = s;
  check_econd11(s, <= 0.0, mcerr);
}

Member Function Documentation

get_integ_active()

double Heed::PointsRan::get_integ_active ( ) const

inline

Definition at line 68 of file PointsRan.h.

{ return integ_active; }

get_integ_total()

double Heed::PointsRan::get_integ_total ( ) const

inline

Definition at line 66 of file PointsRan.h.

{ return integ_total; }

print()

void Heed::PointsRan::print

(

std::ostream &

file

)

const

Definition at line 120 of file PointsRan.cpp.

                                            {
  Ifile << "PointsRan:\n";
  indn.n += 2;
  Ifile << "xmin=" << xmin << " xmax=" << xmax << '\n';
  Ifile << "n_start=" << n_start << " n_finish=" << n_finish << '\n';
  Ifile << "integ_start=" << integ_start << " integ_finish=" << integ_finish
        << '\n';
  Ifile << "integ_total=" << integ_total << " integ_active=" << integ_active
        << '\n';
  // Iprintn(file, integ);
  const long q = x.size();
  Iprintn(file, q);
  for (long n = 0; n < q; n++) {
    file << std::setw(3) << n << ' ' << std::setw(12) << x[n] << ' '
         << std::setw(12) << y[n] << ' ' << std::setw(12) << iy[n];
    if (n < q - 1) file << ' ' << std::setw(12) << a[n];
    file << '\n';
  }
  indn.n -= 2;
}

Referenced by Heed::PairProd::print().

ran()

double Heed::PointsRan::ran

(

double

flat_ran

)

const

Definition at line 91 of file PointsRan.cpp.

                                           {
  mfunnamep("double PointsRan::ran(double flat_ran) const");
  flat_ran = integ_start + integ_active * flat_ran;
  // long q = x.get_qel();
  long n1 = n_start;
  long n2 = n_finish;
  long n3;
  while (n2 - n1 > 1) {
    n3 = n1 + (n2 - n1) / 2;
    if (flat_ran < iy[n3]) {
      n2 = n3;
    } else {
      n1 = n3;
    }
  }
  double dran = flat_ran - iy[n1];
  // double dx = sqrt(2.0 * dran);
  double dx;
  if (a[n1] != 0.0) {
    dx = (-y[n1] + sqrt(y[n1] * y[n1] + 2.0 * a[n1] * dran)) / a[n1];
  } else {
    dx = (x[n2] - x[n1]) / (iy[n2] - iy[n1]) * dran;
  }
  // check_econd11(dx , < 0 , mcerr); // for debug
  // check_econd11(dx , > x[n2] - x[n1] , mcerr); // for debug
  double r = x[n1] + dx;
  return r;
}

Referenced by Heed::PairProd::get_eloss().

---

The documentation for this class was generated from the following files:

• PointsRan.h
• PointsRan.cpp