PointsRan Class Reference

#include <PointsRan.h>

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

Detailed Description

Definition at line 25 of file PointsRan.h.

Constructor & Destructor Documentation

PointsRan() [1/2]

PointsRan::PointsRan ( void  )

inline

Definition at line 43 of file PointsRan.h.

{ ; }

PointsRan() [2/2]

PointsRan::PointsRan	(	DynLinArr< double >	fx,
		DynLinArr< double >	fy,
		double	fxmin,
		double	fxmax
	)

Definition at line 16 of file PointsRan.cpp.

    : xmin(fxmin), xmax(fxmax), x(fx), y(fy) {
  mfunnamep("PointsRan::PointsRan(...)");
  check_econd12(x.get_qel(), !=, y.get_qel(), mcerr);
  check_econd11(x.get_qel(), < 2, mcerr);
  check_econd12(xmin, >=, xmax, mcerr);
  //check_econd12(xmin , > , x[0] , mcerr);
  long q = x.get_qel();
  //check_econd12(xmax , < , x[q-1] , mcerr);
  long n;
  for (n = 0; n < q - 1; n++) {
    check_econd12(x[n], >=, x[n + 1], mcerr);
  }
  for (n = 0; n < q; n++) {
    check_econd11(y[n], < 0.0, mcerr);
  }
  iy = DynLinArr<double>(q);
  a = DynLinArr<double>(q - 1);
  for (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 (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;
  /*
  for( n=0; n<q; n++)
  {
    mcout<<std::setw(3)<<n
        <<' '<<std::setw(12)<<x[n]
        <<' '<<std::setw(12)<<y[n]
        <<' '<<std::setw(12)<<iy[n];
    if(n < q-1)
      mcout<<' '<<std::setw(12)<<a[n];
    mcout<<'\n';
  }
  */
  double s = iy[q - 1];
  integ_total = s;
  check_econd11(s, <= 0.0, mcerr);
 
  //s = 1.0 / s;
  //for( n=0; n<q; n++)
  //{
  //  y[n] = y[n] * s;
  //  iy[n] = iy[n] * s;
  //  if(n < q-1)
  //    a[n] = a[n] * s;
  //}
  //check_econd11(fabs(iy[q-1] - 1.0) , > 1.0e-10 , mcerr);
  //iy[q-1] = 1.0;
}

Member Function Documentation

get_integ_active()

double PointsRan::get_integ_active ( void  ) const

inline

Definition at line 57 of file PointsRan.h.

{ return integ_active; }

get_integ_total()

double PointsRan::get_integ_total ( void  ) const

inline

Definition at line 56 of file PointsRan.h.

{ return integ_total; }

print()

void PointsRan::print

(

std::ostream &

file

)

const

Definition at line 143 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);
  long q = x.get_qel();
  Iprintn(file, q);
  long n;
  for (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 PointsRan::ran

(

double

flat_ran

)

const

Definition at line 114 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::ElElasticScat::fill_hist_low_scat(), and Heed::PairProd::get_eloss().

---

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

• PointsRan.h
• PointsRan.cpp