Geant4_10
RandGaussZiggurat.h
Go to the documentation of this file.
1 /*
2 Code adapted from:
3 http://www.jstatsoft.org/v05/i08/
4 
5 
6 Original disclaimer:
7 
8 The ziggurat method for RNOR and REXP
9 Combine the code below with the main program in which you want
10 normal or exponential variates. Then use of RNOR in any expression
11 will provide a standard normal variate with mean zero, variance 1,
12 while use of REXP in any expression will provide an exponential variate
13 with density exp(-x),x>0.
14 Before using RNOR or REXP in your main, insert a command such as
15 zigset(86947731 );
16 with your own choice of seed value>0, rather than 86947731.
17 (If you do not invoke zigset(...) you will get all zeros for RNOR and REXP.)
18 For details of the method, see Marsaglia and Tsang, "The ziggurat method
19 for generating random variables", Journ. Statistical Software.
20 
21 */
22 
23 
24 #ifndef RandGaussZiggurat_h
25 #define RandGaussZiggurat_h 1
26 
27 #include "cmath"
28 #include "CLHEP/Random/RandGauss.h"
29 
30 namespace CLHEP {
31 
36 class RandGaussZiggurat : public RandGauss {
37 
38 public:
39 
40  inline RandGaussZiggurat ( HepRandomEngine& anEngine, double mean=0.0, double stdDev=1.0 );
41  inline RandGaussZiggurat ( HepRandomEngine* anEngine, double mean=0.0, double stdDev=1.0 );
42 
43  // Destructor
44  virtual ~RandGaussZiggurat();
45 
46  // Static methods to shoot random values using the static generator
47 
48  static inline float shoot() {return ziggurat_RNOR(HepRandom::getTheEngine());};
49  static inline float shoot( float mean, float stdDev ) {return shoot()*stdDev + mean;};
50 
51  static void shootArray ( const int size, float* vect, float mean=0.0, float stdDev=1.0 );
52  static void shootArray ( const int size, double* vect, double mean=0.0, double stdDev=1.0 );
53 
54  // Static methods to shoot random values using a given engine
55  // by-passing the static generator.
56 
57  static inline float shoot( HepRandomEngine* anotherEngine ) {return ziggurat_RNOR(anotherEngine);};
58  static inline float shoot( HepRandomEngine* anotherEngine, float mean, float stdDev ) {return shoot(anotherEngine)*stdDev + mean;};
59 
60  static void shootArray ( HepRandomEngine* anotherEngine, const int size, float* vect, float mean=0.0, float stdDev=1.0 );
61  static void shootArray ( HepRandomEngine* anotherEngine, const int size, double* vect, double mean=0.0, double stdDev=1.0 );
62 
63  // Instance methods using the localEngine to instead of the static
64  // generator, and the default mean and stdDev established at construction
65 
66  inline float fire() {return ziggurat_RNOR(localEngine.get()) * defaultStdDev + defaultMean;};
67 
68  inline float fire ( float mean, float stdDev ) {return ziggurat_RNOR(localEngine.get()) * stdDev + mean;};
69 
70  void fireArray ( const int size, float* vect);
71  void fireArray ( const int size, double* vect);
72  void fireArray ( const int size, float* vect, float mean, float stdDev );
73  void fireArray ( const int size, double* vect, double mean, double stdDev );
74 
75  virtual double operator()();
76  virtual double operator()( double mean, double stdDev );
77 
78  // Save and restore to/from streams
79 
80  std::ostream & put ( std::ostream & os ) const;
81  std::istream & get ( std::istream & is );
82 
83  std::string name() const;
85 
86  static std::string distributionName() {return "RandGaussZiggurat";}
87  // Provides the name of this distribution class
88 
89  static bool ziggurat_init();
90 protected:
92  // Ziggurat Original code:
94  //static unsigned long jz,jsr=123456789;
95  //
96  //#define SHR3 (jz=jsr, jsr^=(jsr<<13), jsr^=(jsr>>17), jsr^=(jsr<<5),jz+jsr)
97  //#define UNI (.5 + (signed) SHR3*.2328306e-9)
98  //#define IUNI SHR3
99  //
100  //static long hz;
101  //static unsigned long iz, kn[128], ke[256];
102  //static float wn[128],fn[128], we[256],fe[256];
103  //
104  //#define RNOR (hz=SHR3, iz=hz&127, (fabs(hz)<kn[iz])? hz*wn[iz] : nfix())
105  //#define REXP (jz=SHR3, iz=jz&255, ( jz <ke[iz])? jz*we[iz] : efix())
106 
107  static unsigned long kn[128], ke[256];
108  static float wn[128],fn[128], we[256],fe[256];
109 
110  static bool ziggurat_is_init;
111 
112  static inline unsigned long ziggurat_SHR3(HepRandomEngine* anEngine) {return (unsigned int)(*anEngine);};
113  static inline float ziggurat_UNI(HepRandomEngine* anEngine) {return anEngine->flat();};
114  static inline float ziggurat_RNOR(HepRandomEngine* anEngine) {
115  long hz=(signed)ziggurat_SHR3(anEngine);
116  unsigned long iz=hz&127;
117  return ((unsigned long)abs(hz)<kn[iz]) ? hz*wn[iz] : ziggurat_nfix(hz,anEngine);
118  };
119  static float ziggurat_nfix(long hz,HepRandomEngine* anEngine);
120 
121 private:
122 
123  // Private copy constructor. Defining it here disallows use.
125 
126  // All the engine info, and the default mean and sigma, are in the RandGauss
127  // base class.
128 
129 };
130 
131 } // namespace CLHEP
132 
133 namespace CLHEP {
134 
135 RandGaussZiggurat::RandGaussZiggurat(HepRandomEngine & anEngine, double mean,double stdDev ): RandGauss(anEngine, mean, stdDev)
136 {
138 }
139 
140 RandGaussZiggurat::RandGaussZiggurat(HepRandomEngine * anEngine, double mean,double stdDev ): RandGauss(anEngine, mean, stdDev)
141 {
143 }
144 
145 } // namespace CLHEP
146 
147 #endif
Float_t d
Definition: plot.C:237
static float ziggurat_nfix(long hz, HepRandomEngine *anEngine)
double defaultStdDev
Definition: RandGauss.h:152
float fire(float mean, float stdDev)
static unsigned long ke[256]
static float shoot(HepRandomEngine *anotherEngine)
virtual double flat()=0
static float ziggurat_UNI(HepRandomEngine *anEngine)
RandGaussZiggurat(HepRandomEngine &anEngine, double mean=0.0, double stdDev=1.0)
static std::string distributionName()
static unsigned long ziggurat_SHR3(HepRandomEngine *anEngine)
static HepRandomEngine * getTheEngine()
Definition: Random.cc:165
static float shoot(float mean, float stdDev)
shared_ptr< HepRandomEngine > localEngine
Definition: RandGauss.h:154
G4double iz
Definition: TRTMaterials.hh:39
void fireArray(const int size, float *vect)
double defaultMean
Definition: RandGauss.h:151
static float ziggurat_RNOR(HepRandomEngine *anEngine)
std::string name() const
std::ostream & put(std::ostream &os) const
static float shoot(HepRandomEngine *anotherEngine, float mean, float stdDev)
static void shootArray(const int size, float *vect, float mean=0.0, float stdDev=1.0)
static unsigned long kn[128]
HepRandomEngine & engine()