Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RandGaussZiggurat.cc
Go to the documentation of this file.
3 #include <iostream>
4 #include <cmath> // for std::log()
5 
6 namespace CLHEP {
7 
11 
13 
15 }
16 
17 std::string RandGaussZiggurat::name() const
18 {
19  return "RandGaussZiggurat";
20 }
21 
23 {
24  const double rzm1 = 2147483648.0, rzm2 = 4294967296.;
25  double dn=3.442619855899,tn=dn,vn=9.91256303526217e-3, q;
26  double de=7.697117470131487, te=de, ve=3.949659822581572e-3;
27  int i;
28 
29 /* Set up tables for RNOR */
30  q=vn/std::exp(-.5*dn*dn);
31  kn[0]=(unsigned long)((dn/q)*rzm1);
32  kn[1]=0;
33 
34  wn[0]=q/rzm1;
35  wn[127]=dn/rzm1;
36 
37  fn[0]=1.;
38  fn[127]=std::exp(-.5*dn*dn);
39 
40  for(i=126;i>=1;i--) {
41  dn=std::sqrt(-2.*std::log(vn/dn+std::exp(-.5*dn*dn)));
42  kn[i+1]=(unsigned long)((dn/tn)*rzm1);
43  tn=dn;
44  fn[i]=std::exp(-.5*dn*dn);
45  wn[i]=dn/rzm1;
46  }
47 
48 /* Set up tables for REXP */
49  q = ve/std::exp(-de);
50  ke[0]=(unsigned long)((de/q)*rzm2);
51  ke[1]=0;
52 
53  we[0]=q/rzm2;
54  we[255]=de/rzm2;
55 
56  fe[0]=1.;
57  fe[255]=std::exp(-de);
58 
59  for(i=254;i>=1;i--) {
60  de=-std::log(ve/de+std::exp(-de));
61  ke[i+1]= (unsigned long)((de/te)*rzm2);
62  te=de;
63  fe[i]=std::exp(-de);
64  we[i]=de/rzm2;
65  }
66  ziggurat_is_init=true;
67 
68  //std::cout<<"Done RandGaussZiggurat::ziggurat_init()"<<std::endl;
69 
70  return true;
71 }
72 
74 {
76  const float r = 3.442620f; /* The start of the right tail */
77  float x, y;
78  unsigned long iz=hz&127;
79  for(;;)
80  {
81  x=hz*wn[iz]; /* iz==0, handles the base strip */
82  if(iz==0) {
83  do {
84  /* change to (1.0 - UNI) as argument to std::log(), because CLHEP generates [0,1),
85  while the original UNI generates (0,1] */
86  x=-std::log(1.0 - ziggurat_UNI(anEngine))*0.2904764; /* .2904764 is 1/r */
87  y=-std::log(1.0 - ziggurat_UNI(anEngine));
88  } while(y+y<x*x);
89  return (hz>0)? r+x : -r-x;
90  }
91  /* iz>0, handle the wedges of other strips */
92  if( fn[iz]+(1.0 - ziggurat_UNI(anEngine))*(fn[iz-1]-fn[iz]) < std::exp(-.5*x*x) ) return x;
93 
94  /* initiate, try to exit for(;;) for loop*/
95  hz=(signed)ziggurat_SHR3(anEngine);
96  iz=hz&127;
97  if((unsigned long)std::abs(hz)<kn[iz]) return (hz*wn[iz]);
98  }
99 }
100 
103 }
104 
105 double RandGaussZiggurat::operator()( double mean, double stdDev ) {
106  return ziggurat_RNOR(localEngine.get()) * stdDev + mean;
107 }
108 
109 void RandGaussZiggurat::shootArray( const int size, float* vect, float mean, float stdDev )
110 {
111  for (int i=0; i<size; ++i) {
112  vect[i] = shoot(mean,stdDev);
113  }
114 }
115 
116 void RandGaussZiggurat::shootArray( const int size, double* vect, double mean, double stdDev )
117 {
118  for (int i=0; i<size; ++i) {
119  vect[i] = shoot(mean,stdDev);
120  }
121 }
122 
123 void RandGaussZiggurat::shootArray( HepRandomEngine* anEngine, const int size, float* vect, float mean, float stdDev )
124 {
125  for (int i=0; i<size; ++i) {
126  vect[i] = shoot(anEngine,mean,stdDev);
127  }
128 }
129 
130 void RandGaussZiggurat::shootArray( HepRandomEngine* anEngine, const int size, double* vect, double mean, double stdDev )
131 {
132  for (int i=0; i<size; ++i) {
133  vect[i] = shoot(anEngine,mean,stdDev);
134  }
135 }
136 
137 void RandGaussZiggurat::fireArray( const int size, float* vect)
138 {
139  for (int i=0; i<size; ++i) {
140  vect[i] = fire( defaultMean, defaultStdDev );
141  }
142 }
143 
144 void RandGaussZiggurat::fireArray( const int size, double* vect)
145 {
146  for (int i=0; i<size; ++i) {
147  vect[i] = fire( defaultMean, defaultStdDev );
148  }
149 }
150 
151 void RandGaussZiggurat::fireArray( const int size, float* vect, float mean, float stdDev )
152 {
153  for (int i=0; i<size; ++i) {
154  vect[i] = fire( mean, stdDev );
155  }
156 }
157 
158 void RandGaussZiggurat::fireArray( const int size, double* vect, double mean, double stdDev )
159 {
160  for (int i=0; i<size; ++i) {
161  vect[i] = fire( mean, stdDev );
162  }
163 }
164 
165 std::ostream & RandGaussZiggurat::put ( std::ostream & os ) const {
166  int pr=os.precision(20);
167  os << " " << name() << "\n";
168  RandGauss::put(os);
169  os.precision(pr);
170  return os;
171 }
172 
173 std::istream & RandGaussZiggurat::get ( std::istream & is ) {
174  std::string inName;
175  is >> inName;
176  if (inName != name()) {
177  is.clear(std::ios::badbit | is.rdstate());
178  std::cerr << "Mismatch when expecting to read state of a "
179  << name() << " distribution\n"
180  << "Name found was " << inName
181  << "\nistream is left in the badbit state\n";
182  return is;
183  }
184  RandGauss::get(is);
185  return is;
186 }
187 
188 } // namespace CLHEP
static CLHEP_THREAD_LOCAL float fe[256]
static CLHEP_THREAD_LOCAL float wn[128]
static float ziggurat_nfix(long hz, HepRandomEngine *anEngine)
double defaultStdDev
Definition: RandGauss.h:153
HepRandomEngine & engine()
Definition: RandGauss.cc:43
tuple x
Definition: test.py:50
static float ziggurat_UNI(HepRandomEngine *anEngine)
std::shared_ptr< HepRandomEngine > localEngine
Definition: RandGauss.h:155
static unsigned long ziggurat_SHR3(HepRandomEngine *anEngine)
std::istream & get(std::istream &is)
std::istream & get(std::istream &is)
Definition: RandGauss.cc:275
static CLHEP_THREAD_LOCAL unsigned long kn[128]
void fireArray(const int size, float *vect)
double defaultMean
Definition: RandGauss.h:152
#define CLHEP_THREAD_LOCAL
Definition: thread_local.h:25
static CLHEP_THREAD_LOCAL unsigned long ke[256]
static float ziggurat_RNOR(HepRandomEngine *anEngine)
std::string name() const
std::ostream & put(std::ostream &os) const
std::ostream & put(std::ostream &os) const
Definition: RandGauss.cc:256
static void shootArray(const int size, float *vect, float mean=0.0, float stdDev=1.0)
static CLHEP_THREAD_LOCAL float fn[128]
static CLHEP_THREAD_LOCAL bool ziggurat_is_init
static CLHEP_THREAD_LOCAL float we[256]
HepRandomEngine & engine()