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 
357 
361  RndVector::iterator nextNumber;
362 
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 */