00001
00002
00003
00004
00005
00006
00007
00008
00009 #ifndef ThePEG_RandomGenerator_H
00010 #define ThePEG_RandomGenerator_H
00011
00012
00013 #include "ThePEG/Config/ThePEG.h"
00014 #include "ThePEG/Interface/Interfaced.h"
00015 #include "gsl/gsl_rng.h"
00016
00017 namespace ThePEG {
00018
00037 class RandomGenerator: public Interfaced {
00038
00039 public:
00040
00042 typedef vector<double> RndVector;
00043
00045 typedef RndVector::size_type size_type;
00046
00047 public:
00048
00054 RandomGenerator();
00055
00059 RandomGenerator(const RandomGenerator &);
00060
00064 virtual ~RandomGenerator();
00066
00071 virtual void setSeed(long seed) = 0;
00072
00079 double rnd() {
00080 if ( nextNumber == theNumbers.end() ) fill();
00081 return *nextNumber++;
00082 }
00083
00088 template <typename Unit> Unit rnd(Unit b) { return b*rnd(); }
00089
00094 template <typename Unit>
00095 Unit rnd(Unit a, Unit b) { return a + rnd(b - a); }
00096
00101 RndVector rndvec(int n) {
00102 RndVector ret(n);
00103 for ( int i = 0; i < n; ++i ) ret[i] = rnd();
00104 return ret;
00105 }
00106
00111 double operator()() { return rnd(); }
00112
00117 long operator()(long N) { return long(rnd() * N); }
00118
00123 bool rndbool(double p = 0.5);
00124
00129 bool rndbool(double p1, double p2) { return rndbool(p1/(p1 + p2)); }
00130
00135 int rndsign(double p1, double p2, double p3);
00136
00141 int rnd2(double p0, double p1) {
00142 return rndbool(p0, p1)? 0: 1;
00143 }
00144
00149 int rnd3(double p0, double p1, double p2) {
00150 return 1 + rndsign(p0, p1, p2);
00151 }
00152
00157 int rnd4(double p0, double p1, double p2, double p3);
00158
00163 double rndExp() {
00164 return -log(rnd());
00165 }
00166
00171 template <typename Unit>
00172 Unit rndExp(Unit mean) { return mean*rndExp(); }
00173
00178 double rndGauss() {
00179 if ( gaussSaved ) {
00180 gaussSaved = false;
00181 return savedGauss;
00182 }
00183 double r = sqrt(-2.0*log(rnd()));
00184 double phi = rnd()*2.0*Constants::pi;
00185 savedGauss = r*cos(phi);
00186 gaussSaved = true;
00187 return r*sin(phi);
00188 }
00189
00194 template <typename Unit>
00195 Unit rndGauss(Unit sigma, Unit mean = Unit()) {
00196 return mean + sigma*rndGauss();
00197 }
00198
00204 template <typename Unit>
00205 Unit rndBW(Unit mean, Unit gamma) {
00206 if ( gamma <= Unit() ) return mean;
00207 return mean + 0.5*gamma*tan(rnd(atan(-2.0*mean/gamma), Constants::pi/2));
00208 }
00209
00216 template <typename Unit>
00217 Unit rndBW(Unit mean, Unit gamma, Unit cut) {
00218 if ( gamma <= Unit() || cut <= Unit() ) return mean;
00219 return mean + 0.5*gamma*tan(rnd(atan(-2.0*min(mean,cut)/gamma),
00220 atan(2.0*cut/gamma)));
00221 }
00222
00227 template <typename Unit>
00228 Unit rndRelBW(Unit mean, Unit gamma) {
00229 if ( gamma <= Unit() ) return mean;
00230 return sqrt(sqr(mean) + mean*gamma*tan(rnd(atan(-mean/gamma),
00231 Constants::pi/2)));
00232 }
00233
00240 template <typename Unit>
00241 Unit rndRelBW(Unit mean, Unit gamma, Unit cut) {
00242 if ( gamma <= Unit() || cut <= Unit() ) return mean;
00243 double minarg = cut > mean? -mean/gamma:
00244 (sqr(mean - cut) - sqr(mean))/(gamma*mean);
00245 double maxarg = (sqr(mean + cut) - sqr(mean))/(mean*gamma);
00246 return sqrt(sqr(mean) + mean*gamma*tan(rnd(atan(minarg), atan(maxarg))));
00247 }
00248
00256 long rndPoisson(double mean);
00258
00275 void push_back(double r) {
00276 if ( r > 0.0 && r < 1.0 && nextNumber != theNumbers.begin() )
00277 *--nextNumber = r;
00278 }
00279
00283 void pop_back() {
00284 if ( nextNumber != theNumbers.end() ) ++nextNumber;
00285 }
00286
00291 void flush() { nextNumber = theNumbers.end(); }
00292
00297 template <typename OutputIterator>
00298 void rnd(OutputIterator o, size_type n) {
00299 while ( n-- ) *o++ = rnd();
00300 }
00302
00303 protected:
00304
00310 virtual void doinit();
00311
00312 public:
00313
00314
00321 void persistentOutput(PersistentOStream & os) const;
00322
00328 void persistentInput(PersistentIStream & is, int version);
00330
00334 static void Init();
00335
00339 gsl_rng * getGslInterface() { return gsl; }
00340
00341 protected:
00342
00346 void setSize(size_type newSize);
00347
00351 virtual void fill() = 0;
00352
00356 RndVector theNumbers;
00357
00361 RndVector::iterator nextNumber;
00362
00366 size_type theSize;
00367
00371 long theSeed;
00372
00376 mutable double savedGauss;
00377
00381 mutable bool gaussSaved;
00382
00386 gsl_rng * gsl;
00387
00388 private:
00389
00394 static ClassDescription<RandomGenerator> initRandomGenerator;
00395
00399 RandomGenerator & operator=(const RandomGenerator &);
00400
00401 };
00402
00407 template <>
00408 struct BaseClassTrait<RandomGenerator,1>: public ClassTraitsType {
00410 typedef Interfaced NthBase;
00411 };
00412
00415 template <>
00416 struct ClassTraits<RandomGenerator>:
00417 public ClassTraitsBase<RandomGenerator> {
00419 static string className() { return "ThePEG::RandomGenerator";
00420 }
00423 static TPtr create() {
00424 throw std::logic_error("Tried to instantiate abstract class"
00425 "'Pythis7::RandomGenerator'");
00426 }
00427 };
00428
00431 }
00432
00433 #endif