RandomGenerator.h

// -*- C++ -*-
//
// RandomGenerator.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2019 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_RandomGenerator_H
#define ThePEG_RandomGenerator_H
// This is the declaration of the RandomGenerator class.

#include "ThePEG/Config/ThePEG.h"
#include "ThePEG/Interface/Interfaced.h"
#include "gsl/gsl_rng.h"

namespace ThePEG {

class RandomGenerator: public Interfaced {

public:

  typedef vector<double> RndVector;

  typedef RndVector::size_type size_type;

public:

  RandomGenerator();

  RandomGenerator(const RandomGenerator &);

  virtual ~RandomGenerator();

  virtual void setSeed(long seed) = 0;

  double rnd() {
    if ( nextNumber == theNumbers.end() ) fill();
    return *nextNumber++;
  }

  template <typename Unit> Unit rnd(Unit b) { return b*rnd(); }

  template <typename Unit>
  Unit rnd(Unit a, Unit b) { return a + rnd(b - a); }

  RndVector rndvec(int n) {
    RndVector ret(n);
    for ( int i = 0; i < n; ++i ) ret[i] = rnd();
    return ret;
  }

  double operator()() { return rnd(); }

  double operator()(double N) { return double(rnd() * N); }
  
  long operator()(long N) { return long(rnd() * N); }

  bool rndbool(double p = 0.5) {
    return rnd() < p;
  }

  bool prndbool(double p = 0.5);

  bool rndbool(double p1, double p2) { return rndbool(p1/(p1 + p2)); }

  bool prndbool(double p1, double p2) { return prndbool(p1/(p1 + p2)); }

  int rndsign(double p1, double p2, double p3);

  int prndsign(double p1, double p2, double p3);

  int rnd2(double p0, double p1) {
    return rndbool(p0, p1)? 0: 1;
  }

  int rnd3(double p0, double p1, double p2) {
    return 1 + rndsign(p0, p1, p2);
  }

  int rnd4(double p0, double p1, double p2, double p3);

  int rnd5(double p0, double p1, double p2, double p3, double p4);  

  double rndExp() {
    return -log(rnd());
  }

  template <typename Unit>
  Unit rndExp(Unit mean) { return mean*rndExp(); }

  void rndGaussTwoNumbers(double & randomNumberOne, double & randomNumberTwo) { 
    double r = sqrt(-2.0*log(rnd()));
    double phi = rnd()*2.0*Constants::pi;
    randomNumberOne = r*sin(phi);
    randomNumberTwo = r*cos(phi);
  }
  
  double rndGauss() { 
    if ( gaussSaved ) {
      gaussSaved = false;
      return savedGauss;
    }
    double  randomNumberOne, randomNumberTwo;
    rndGaussTwoNumbers(randomNumberOne, randomNumberTwo);
    savedGauss = randomNumberTwo;
    gaussSaved = true;
    return randomNumberOne;
  }

  template <typename Unit>
  Unit rndGauss(Unit sigma, Unit mean = Unit()) {
    return mean + sigma*rndGauss();
  }

  template <typename Unit>
  void rndGaussTwoNumbers(Unit & randomNumberOne, Unit & randomNumberTwo, Unit sigma, Unit mean = Unit()) {
    double r1,r2;
    rndGaussTwoNumbers(r1,r2);
    randomNumberOne = mean + sigma * r1;
    randomNumberTwo = mean + sigma * r2;
  }
  
  template <typename Unit>
  Unit rndBW(Unit mean, Unit gamma) {
    if ( gamma <= Unit() ) return mean;
    return mean + 0.5*gamma*tan(rnd(atan(-2.0*mean/gamma), Constants::pi/2));
  }

  template <typename Unit>
  Unit rndBW(Unit mean, Unit gamma, Unit cut) {
    if ( gamma <= Unit() || cut <= Unit() ) return mean;
    return mean + 0.5*gamma*tan(rnd(atan(-2.0*min(mean,cut)/gamma),
                                    atan(2.0*cut/gamma)));
  }

  template <typename Unit>
  Unit rndRelBW(Unit mean, Unit gamma) {
    if ( gamma <= Unit() ) return mean;
    return sqrt(sqr(mean) + mean*gamma*tan(rnd(atan(-mean/gamma),
                                               Constants::pi/2)));
  }

  template <typename Unit>
  Unit rndRelBW(Unit mean, Unit gamma, Unit cut) {
    if ( gamma <= Unit() || cut <= Unit() ) return mean;
    double minarg = cut > mean? -mean/gamma:
      (sqr(mean - cut) - sqr(mean))/(gamma*mean);
    double maxarg = (sqr(mean + cut) - sqr(mean))/(mean*gamma);
    return sqrt(sqr(mean) + mean*gamma*tan(rnd(atan(minarg), atan(maxarg))));
  }

  long rndPoisson(double mean);

  void push_back(double r) {  
    if ( r > 0.0 && r < 1.0 && nextNumber != theNumbers.begin() )
      *--nextNumber = r;
  }

  void pop_back() {
    if ( nextNumber != theNumbers.end() ) ++nextNumber;
  }

  void flush() {
    nextNumber = theNumbers.end();
    gaussSaved = false;
  }

  template <typename OutputIterator>
  void rnd(OutputIterator o, size_type n) {
    while ( n-- ) *o++ = rnd();
  }

protected:

  virtual void doinit();

public:


  void persistentOutput(PersistentOStream & os) const;

  void persistentInput(PersistentIStream & is, int version);

  static void Init();

  gsl_rng * getGslInterface() { return gsl; }

protected:

  void setSize(size_type newSize);

  virtual void fill() = 0;

  RndVector theNumbers;

  RndVector::iterator nextNumber;

  size_type theSize;

  long theSeed;

  mutable double savedGauss;

  mutable bool gaussSaved;

  gsl_rng * gsl;

private:

  static ClassDescription<RandomGenerator> initRandomGenerator;

  RandomGenerator & operator=(const RandomGenerator &) = delete;

};

template <>
struct BaseClassTrait<RandomGenerator,1>: public ClassTraitsType {
  typedef Interfaced NthBase;
};

template <>
struct ClassTraits<RandomGenerator>:
    public ClassTraitsBase<RandomGenerator> {
  static string className() { return "ThePEG::RandomGenerator";
  }
  static TPtr create() {
    throw std::logic_error("Tried to instantiate abstract class"
                           "'Pythis7::RandomGenerator'");
  }
};

}

#endif /* ThePEG_RandomGenerator_H */