// RandomGenerator.cpp

#include "RandomGenerator.h"
#include <iostream>
#include "CLHEP/Random/RandGauss.h"
#include "trfutil/trfstream.h"

using std::cout;

//**********************************************************************

// Default constructor.

RandomGenerator::RandomGenerator() : _rand() {
  set_seed(19970729);
}

//**********************************************************************

// Constructor from seed.

RandomGenerator::RandomGenerator(long seed) : _rand() {
  set_seed(seed);
}

//**********************************************************************

// Copy constructor.
// Set the original seed and skip through the same number of
// values to obtain the same state.

RandomGenerator::RandomGenerator(const RandomGenerator& rhs)
: _rand() {
  set_state(rhs);
}

//**********************************************************************

// Destructor.

RandomGenerator::~RandomGenerator() { }

//**********************************************************************

// Set seed.
// This reinitializes the generator engine, positioning it at the
// the first value in the sequence.

void RandomGenerator::set_seed(long seed) {
  _ncall = 0;
  _rand.setSeed(seed,3);
}

//**********************************************************************

// Set the state.

void RandomGenerator::set_state(long seed, int ncall) {
  set_seed(seed);
  for ( int i=0; i<ncall; ++i ) flat(0.0,1.0);
}

//**********************************************************************

// Set the state (seed + ncall) to be the same as that of the argument.

void RandomGenerator::set_state(const RandomGenerator& gen) {
  set_state(gen._rand.getSeed(),gen._ncall);
}


//**********************************************************************

// Return a value from a flat distribution in range (min,max).

double RandomGenerator::flat(double min, double max) {
  ++_ncall;
  return min + (max-min)*_rand.flat();
}

//**********************************************************************

// Return a value from a Gaussian distribution with sigma = 1.0.

double RandomGenerator::gauss() {
  double v1, v2, r;
  do {
    v1 = flat( -1.0, 1.0 );
    v2 = flat( -1.0, 1.0 );
    r = v1*v1 + v2*v2;
  } while ( r > 1.0 );
  double fac = sqrt( -2.0*log(r)/r);
  return v2*fac;
}

//**********************************************************************

// output stream

void RandomGenerator::ostr(ostream& stream) const {
  cout << begin_object;
  raw_ostr(stream);
  cout << end_object;
}

//**********************************************************************

// raw output stream

void RandomGenerator:: raw_ostr(ostream& stream) const {
  stream << "Seed: " << _rand.getSeed() << new_line;
  stream << "# calls: " << _ncall;
}

//**********************************************************************

// output stream function

ostream& operator<<(ostream& stream, const RandomGenerator& rgen) {
  rgen.ostr(stream);
  return stream;
}

//**********************************************************************