thepeg is hosted by Hepforge, IPPP Durham
ThePEG  2.1.5
RandomGenerator.h
1 // -*- C++ -*-
2 //
3 // RandomGenerator.h is a part of ThePEG - Toolkit for HEP Event Generation
4 // Copyright (C) 1999-2017 Leif Lonnblad
5 //
6 // ThePEG is licenced under version 3 of the GPL, see COPYING for details.
7 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
8 //
9 #ifndef ThePEG_RandomGenerator_H
10 #define ThePEG_RandomGenerator_H
11 // This is the declaration of the RandomGenerator class.
12 
13 #include "ThePEG/Config/ThePEG.h"
14 #include "ThePEG/Interface/Interfaced.h"
15 #include "gsl/gsl_rng.h"
16 
17 namespace ThePEG {
18 
37 class RandomGenerator: public Interfaced {
38 
39 public:
40 
42  typedef vector<double> RndVector;
43 
45  typedef RndVector::size_type size_type;
46 
47 public:
48 
55 
60 
64  virtual ~RandomGenerator();
66 
71  virtual void setSeed(long seed) = 0;
72 
79  double rnd() {
80  if ( nextNumber == theNumbers.end() ) fill();
81  return *nextNumber++;
82  }
83 
88  template <typename Unit> Unit rnd(Unit b) { return b*rnd(); }
89 
94  template <typename Unit>
95  Unit rnd(Unit a, Unit b) { return a + rnd(b - a); }
96 
101  RndVector rndvec(int n) {
102  RndVector ret(n);
103  for ( int i = 0; i < n; ++i ) ret[i] = rnd();
104  return ret;
105  }
106 
111  double operator()() { return rnd(); }
112 
119  double operator()(double N) { return double(rnd() * N); }
120 
125  long operator()(long N) { return long(rnd() * N); }
126 
130  bool rndbool(double p = 0.5) {
131  return rnd() < p;
132  }
133 
138  bool prndbool(double p = 0.5);
139 
144  bool rndbool(double p1, double p2) { return rndbool(p1/(p1 + p2)); }
145 
150  bool prndbool(double p1, double p2) { return prndbool(p1/(p1 + p2)); }
151 
156  int rndsign(double p1, double p2, double p3);
157 
162  int prndsign(double p1, double p2, double p3);
163 
168  int rnd2(double p0, double p1) {
169  return rndbool(p0, p1)? 0: 1;
170  }
171 
176  int rnd3(double p0, double p1, double p2) {
177  return 1 + rndsign(p0, p1, p2);
178  }
179 
184  int rnd4(double p0, double p1, double p2, double p3);
185 
190  double rndExp() {
191  return -log(rnd());
192  }
193 
198  template <typename Unit>
199  Unit rndExp(Unit mean) { return mean*rndExp(); }
200 
208  void rndGaussTwoNumbers(double & randomNumberOne, double & randomNumberTwo) {
209  double r = sqrt(-2.0*log(rnd()));
210  double phi = rnd()*2.0*Constants::pi;
211  randomNumberOne = r*sin(phi);
212  randomNumberTwo = r*cos(phi);
213  }
214 
221  double rndGauss() {
222  if ( gaussSaved ) {
223  gaussSaved = false;
224  return savedGauss;
225  }
226  double randomNumberOne, randomNumberTwo;
227  rndGaussTwoNumbers(randomNumberOne, randomNumberTwo);
228  savedGauss = randomNumberTwo;
229  gaussSaved = true;
230  return randomNumberOne;
231  }
232 
237  template <typename Unit>
238  Unit rndGauss(Unit sigma, Unit mean = Unit()) {
239  return mean + sigma*rndGauss();
240  }
241 
246  template <typename Unit>
247  void rndGaussTwoNumbers(Unit & randomNumberOne, Unit & randomNumberTwo, Unit sigma, Unit mean = Unit()) {
248  double r1,r2;
249  rndGaussTwoNumbers(r1,r2);
250  randomNumberOne = mean + sigma * r1;
251  randomNumberTwo = mean + sigma * r2;
252  }
253 
259  template <typename Unit>
260  Unit rndBW(Unit mean, Unit gamma) {
261  if ( gamma <= Unit() ) return mean;
262  return mean + 0.5*gamma*tan(rnd(atan(-2.0*mean/gamma), Constants::pi/2));
263  }
264 
271  template <typename Unit>
272  Unit rndBW(Unit mean, Unit gamma, Unit cut) {
273  if ( gamma <= Unit() || cut <= Unit() ) return mean;
274  return mean + 0.5*gamma*tan(rnd(atan(-2.0*min(mean,cut)/gamma),
275  atan(2.0*cut/gamma)));
276  }
277 
282  template <typename Unit>
283  Unit rndRelBW(Unit mean, Unit gamma) {
284  if ( gamma <= Unit() ) return mean;
285  return sqrt(sqr(mean) + mean*gamma*tan(rnd(atan(-mean/gamma),
286  Constants::pi/2)));
287  }
288 
295  template <typename Unit>
296  Unit rndRelBW(Unit mean, Unit gamma, Unit cut) {
297  if ( gamma <= Unit() || cut <= Unit() ) return mean;
298  double minarg = cut > mean? -mean/gamma:
299  (sqr(mean - cut) - sqr(mean))/(gamma*mean);
300  double maxarg = (sqr(mean + cut) - sqr(mean))/(mean*gamma);
301  return sqrt(sqr(mean) + mean*gamma*tan(rnd(atan(minarg), atan(maxarg))));
302  }
303 
311  long rndPoisson(double mean);
313 
330  void push_back(double r) {
331  if ( r > 0.0 && r < 1.0 && nextNumber != theNumbers.begin() )
332  *--nextNumber = r;
333  }
334 
338  void pop_back() {
339  if ( nextNumber != theNumbers.end() ) ++nextNumber;
340  }
341 
346  void flush() {
347  nextNumber = theNumbers.end();
348  gaussSaved = false;
349  }
350 
355  template <typename OutputIterator>
356  void rnd(OutputIterator o, size_type n) {
357  while ( n-- ) *o++ = rnd();
358  }
360 
361 protected:
362 
368  virtual void doinit();
369 
370 public:
371 
372 
379  void persistentOutput(PersistentOStream & os) const;
380 
386  void persistentInput(PersistentIStream & is, int version);
388 
392  static void Init();
393 
397  gsl_rng * getGslInterface() { return gsl; }
398 
399 protected:
400 
404  void setSize(size_type newSize);
405 
409  virtual void fill() = 0;
410 
414  RndVector theNumbers;
415 
419  RndVector::iterator nextNumber;
420 
424  size_type theSize;
425 
429  long theSeed;
430 
434  mutable double savedGauss;
435 
439  mutable bool gaussSaved;
440 
444  gsl_rng * gsl;
445 
446 private:
447 
453 
457  RandomGenerator & operator=(const RandomGenerator &) = delete;
458 
459 };
460 
465 template <>
468  typedef Interfaced NthBase;
469 };
470 
473 template <>
474 struct ClassTraits<RandomGenerator>:
475  public ClassTraitsBase<RandomGenerator> {
477  static string className() { return "ThePEG::RandomGenerator";
478  }
481  static TPtr create() {
482  throw std::logic_error("Tried to instantiate abstract class"
483  "'Pythis7::RandomGenerator'");
484  }
485 };
486 
489 }
490 
491 #endif /* ThePEG_RandomGenerator_H */
RndVector rndvec(int n)
Return n flat random number in the interval .
gsl_rng * gsl
A pinter to a gsl_rng interface to this generator.
PersistentIStream is used to read persistent objects from a stream where they were previously written...
Unit rnd(Unit a, Unit b)
Return a flat random number in the interval .
void flush()
Discard all random numbers in the cache.
ClassTraitsType is an empty, non-polymorphic, base class.
Definition: ClassTraits.h:30
A concreate implementation of ClassDescriptionBase describing a concrete class with persistent data...
RndVector::iterator nextNumber
Iterator pointing to the next number to be extracted.
RandomGenerator()
Default constructor.
PersistentOStream is used to write objects persistently to a stream from which they can be read in ag...
size_type theSize
The size of the cache.
void push_back(double r)
Give back a partly unused random number.
long operator()(long N)
Return a (possibly cached) flat integer random number in the interval .
constexpr double pi
Good old .
Definition: Constants.h:54
static void Init()
Standard Init function used to initialize the interface.
double operator()(double N)
Return a (possibly cached) flat integer random number in the interval .
double rndGauss()
Return a number distributed according to a Gaussian distribution with zero mean and unit variance...
virtual void setSeed(long seed)=0
Reset the underlying random engine with the given seed.
This is the main namespace within which all identifiers in ThePEG are declared.
Definition: FactoryBase.h:28
RandomGenerator is an interface to the CLHEP::RandomEngine classes.
This is the main config header file for ThePEG.
Unit rndBW(Unit mean, Unit gamma)
Return a positive number distributed according to a non-relativistic Breit-Wigner with a given width...
Unit rndExp(Unit mean)
Return a number between zero and infinity, distributed according to where is the mean value...
int rnd2(double p0, double p1)
Return an integer with probability p /(p0+p1).
int rnd3(double p0, double p1, double p2)
Return an integer with probability p /(p0+p1+p2).
double rnd()
Return a (possibly cached) flat random number in the interval .
static ClassDescription< RandomGenerator > initRandomGenerator
Describe a concrete class with persistent data.
RandomGenerator & operator=(const RandomGenerator &)=delete
Private and non-existent assignment operator.
void pop_back()
Discard the next random number in the cache.
void setSize(size_type newSize)
Utility function for the interface.
virtual void fill()=0
Fill the cache with random numbers.
Unit rnd(Unit b)
Return a flat random number in the interval .
virtual void doinit()
Initializes this random generator.
RndVector theNumbers
A vector of cached numbers.
double operator()()
Return a (possibly cached) flat random number in the interval .
void persistentInput(PersistentIStream &is, int version)
Function used to read in object persistently.
Unit rndRelBW(Unit mean, Unit gamma)
Return a positive number distributed according to a relativistic Breit-Wigner with a given width...
void persistentOutput(PersistentOStream &os) const
Function used to write out object persistently.
void rndGaussTwoNumbers(Unit &randomNumberOne, Unit &randomNumberTwo, Unit sigma, Unit mean=Unit())
Return two numbers distributed according to a Gaussian distribution with a given standard deviation...
vector< double > RndVector
A vector of doubles.
int prndsign(double p1, double p2, double p3)
Return -1, 0, or 1 with relative probabilities p1, p2, p3.
int rndsign(double p1, double p2, double p3)
Return -1, 0, or 1 with relative probabilities p1, p2, p3.
bool gaussSaved
Indicate the precense of a saved Gaussian random number.
bool rndbool(double p1, double p2)
Return a true with probability p1/(p1+p2).
The Interfaced class is derived from the InterfacedBase class adding a couple of things particular to...
Definition: Interfaced.h:38
bool rndbool(double p=0.5)
Return a true with probability p.
Unit rndGauss(Unit sigma, Unit mean=Unit())
Return a number distributed according to a Gaussian distribution with a given standard deviation...
Unit rndBW(Unit mean, Unit gamma, Unit cut)
Return a positive number distributed according to a non-relativistic Breit-Wigner with a given width...
double savedGauss
A saved Gaussian random number.
void rnd(OutputIterator o, size_type n)
Generate n random numbers between 0 and 1 and put them in the output iterator.
int rnd4(double p0, double p1, double p2, double p3)
Return an integer/ with probability p (p0+p1+p2+p3).
The default concrete implementation of ClassTraitsBase.
Definition: ClassTraits.h:134
long rndPoisson(double mean)
Return a non-negative number generated according to a Poissonian distribution with a given mean...
Unit rndRelBW(Unit mean, Unit gamma, Unit cut)
Return a positive number distributed according to a relativistic Breit-Wigner with a given width...
bool prndbool(double p=0.5)
Return a true with probability p.
long theSeed
The seed to initialize the random generator with.
double rndExp()
Return a number between zero and infinity, distributed according to .
RndVector::size_type size_type
The size_type of RndVector.
BaseClassTraits describes the base classes of the templated class.
Definition: ClassTraits.h:156
void rndGaussTwoNumbers(double &randomNumberOne, double &randomNumberTwo)
Return two numbers distributed according to a Gaussian distribution with zero mean and unit variance...
gsl_rng * getGslInterface()
Return a gsl_rng interface to this random generator.
virtual ~RandomGenerator()
Destructor.
The templated ClassTraitsBase class defines a set of default information about classes used by ThePEG...
Definition: ClassTraits.h:52
bool prndbool(double p1, double p2)
Return a true with probability p1/(p1+p2).