ThePEG  1.8.0
RandomGenerator.h
1 // -*- C++ -*-
2 //
3 // RandomGenerator.h is a part of ThePEG - Toolkit for HEP Event Generation
4 // Copyright (C) 1999-2011 Leif Lonnblad
5 //
6 // ThePEG is licenced under version 2 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 
117  long operator()(long N) { return long(rnd() * N); }
118 
123  bool rndbool(double p = 0.5);
124 
129  bool rndbool(double p1, double p2) { return rndbool(p1/(p1 + p2)); }
130 
135  int rndsign(double p1, double p2, double p3);
136 
141  int rnd2(double p0, double p1) {
142  return rndbool(p0, p1)? 0: 1;
143  }
144 
149  int rnd3(double p0, double p1, double p2) {
150  return 1 + rndsign(p0, p1, p2);
151  }
152 
157  int rnd4(double p0, double p1, double p2, double p3);
158 
163  double rndExp() {
164  return -log(rnd());
165  }
166 
171  template <typename Unit>
172  Unit rndExp(Unit mean) { return mean*rndExp(); }
173 
178  double rndGauss() {
179  if ( gaussSaved ) {
180  gaussSaved = false;
181  return savedGauss;
182  }
183  double r = sqrt(-2.0*log(rnd()));
184  double phi = rnd()*2.0*Constants::pi;
185  savedGauss = r*cos(phi);
186  gaussSaved = true;
187  return r*sin(phi);
188  }
189 
194  template <typename Unit>
195  Unit rndGauss(Unit sigma, Unit mean = Unit()) {
196  return mean + sigma*rndGauss();
197  }
198 
204  template <typename Unit>
205  Unit rndBW(Unit mean, Unit gamma) {
206  if ( gamma <= Unit() ) return mean;
207  return mean + 0.5*gamma*tan(rnd(atan(-2.0*mean/gamma), Constants::pi/2));
208  }
209 
216  template <typename Unit>
217  Unit rndBW(Unit mean, Unit gamma, Unit cut) {
218  if ( gamma <= Unit() || cut <= Unit() ) return mean;
219  return mean + 0.5*gamma*tan(rnd(atan(-2.0*min(mean,cut)/gamma),
220  atan(2.0*cut/gamma)));
221  }
222 
227  template <typename Unit>
228  Unit rndRelBW(Unit mean, Unit gamma) {
229  if ( gamma <= Unit() ) return mean;
230  return sqrt(sqr(mean) + mean*gamma*tan(rnd(atan(-mean/gamma),
231  Constants::pi/2)));
232  }
233 
240  template <typename Unit>
241  Unit rndRelBW(Unit mean, Unit gamma, Unit cut) {
242  if ( gamma <= Unit() || cut <= Unit() ) return mean;
243  double minarg = cut > mean? -mean/gamma:
244  (sqr(mean - cut) - sqr(mean))/(gamma*mean);
245  double maxarg = (sqr(mean + cut) - sqr(mean))/(mean*gamma);
246  return sqrt(sqr(mean) + mean*gamma*tan(rnd(atan(minarg), atan(maxarg))));
247  }
248 
256  long rndPoisson(double mean);
258 
275  void push_back(double r) {
276  if ( r > 0.0 && r < 1.0 && nextNumber != theNumbers.begin() )
277  *--nextNumber = r;
278  }
279 
283  void pop_back() {
284  if ( nextNumber != theNumbers.end() ) ++nextNumber;
285  }
286 
291  void flush() { nextNumber = theNumbers.end(); }
292 
297  template <typename OutputIterator>
298  void rnd(OutputIterator o, size_type n) {
299  while ( n-- ) *o++ = rnd();
300  }
302 
303 protected:
304 
310  virtual void doinit();
311 
312 public:
313 
314 
321  void persistentOutput(PersistentOStream & os) const;
322 
328  void persistentInput(PersistentIStream & is, int version);
330 
334  static void Init();
335 
339  gsl_rng * getGslInterface() { return gsl; }
340 
341 protected:
342 
346  void setSize(size_type newSize);
347 
351  virtual void fill() = 0;
352 
356  RndVector theNumbers;
357 
361  RndVector::iterator nextNumber;
362 
366  size_type theSize;
367 
371  long theSeed;
372 
376  mutable double savedGauss;
377 
381  mutable bool gaussSaved;
382 
386  gsl_rng * gsl;
387 
388 private:
389 
395 
400 
401 };
402 
407 template <>
410  typedef Interfaced NthBase;
411 };
412 
415 template <>
416 struct ClassTraits<RandomGenerator>:
417  public ClassTraitsBase<RandomGenerator> {
419  static string className() { return "ThePEG::RandomGenerator";
420  }
423  static TPtr create() {
424  throw std::logic_error("Tried to instantiate abstract class"
425  "'Pythis7::RandomGenerator'");
426  }
427 };
428 
431 }
432 
433 #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 .
static void Init()
Standard Init function used to initialize the interface.
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.
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.
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...
vector< double > RndVector
A vector of doubles.
void persistentOutput(PersistentOStream &os) const
Function used to write out object persistently.
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.
const double pi
Good old .
Definition: Constants.h:54
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).
RandomGenerator & operator=(const RandomGenerator &)
Private and non-existent assignment operator.
static TPtr create()
Create a T object and return a smart pointer to it.
Definition: ClassTraits.h:60
long rndPoisson(double mean)
Return a non-negative number generated according to a Poissonian distribution with a given mean...
static string className()
Return the name of class T.
Definition: ClassTraits.h:66
int NthBase
The type of the BaseN'th base class (int means there are no further base classes).
Definition: ClassTraits.h:161
Unit rndRelBW(Unit mean, Unit gamma, Unit cut)
Return a positive number distributed according to a relativistic Breit-Wigner with a given width...
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.
ThePEG::Ptr< T >::pointer TPtr
Alias for a reference counted pointer to T .
Definition: ClassTraits.h:54
BaseClassTraits describes the base classes of the templated class.
Definition: ClassTraits.h:156
gsl_rng * getGslInterface()
Return a gsl_rng interface to this random generator.
virtual ~RandomGenerator()
Destructor.