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 */