ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RandExpZiggurat.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RandExpZiggurat.cc
2 
4 
5 #include <iostream>
6 #include <cmath> // for std::log()
7 
8 namespace CLHEP {
9 
13 
14 std::string RandExpZiggurat::name() const {return "RandExpZiggurat";}
15 
17 
19 }
20 
21 RandExpZiggurat::RandExpZiggurat(const RandExpZiggurat& right) : HepRandom(right),defaultMean(right.defaultMean)
22 {
23 }
24 
26 {
27  return fire( defaultMean );
28 }
29 
30 void RandExpZiggurat::shootArray( const int size, float* vect, float mean )
31 {
32  for (int i=0; i<size; ++i) vect[i] = shoot(mean);
33 }
34 
35 void RandExpZiggurat::shootArray( const int size, double* vect, double mean )
36 {
37  for (int i=0; i<size; ++i) vect[i] = shoot(mean);
38 }
39 
40 void RandExpZiggurat::shootArray(HepRandomEngine* anEngine, const int size, float* vect, float mean )
41 {
42  for (int i=0; i<size; ++i) vect[i] = shoot(anEngine, mean);
43 }
44 
45 void RandExpZiggurat::shootArray(HepRandomEngine* anEngine, const int size, double* vect, double mean )
46 {
47  for (int i=0; i<size; ++i) vect[i] = shoot(anEngine, mean);
48 }
49 
50 void RandExpZiggurat::fireArray( const int size, float* vect)
51 {
52  for (int i=0; i<size; ++i) vect[i] = fire( defaultMean );
53 }
54 
55 void RandExpZiggurat::fireArray( const int size, double* vect)
56 {
57  for (int i=0; i<size; ++i) vect[i] = fire( defaultMean );
58 }
59 
60 void RandExpZiggurat::fireArray( const int size, float* vect, float mean )
61 {
62  for (int i=0; i<size; ++i) vect[i] = fire( mean );
63 }
64 
65 void RandExpZiggurat::fireArray( const int size, double* vect, double mean )
66 {
67  for (int i=0; i<size; ++i) vect[i] = fire( mean );
68 }
69 
70 std::ostream & RandExpZiggurat::put ( std::ostream & os ) const {
71  int pr=os.precision(20);
72  std::vector<unsigned long> t(2);
73  os << " " << name() << "\n";
74  os << "Uvec" << "\n";
76  os << defaultMean << " " << t[0] << " " << t[1] << "\n";
77  os.precision(pr);
78  return os;
79 #ifdef REMOVED
80  int pr=os.precision(20);
81  os << " " << name() << "\n";
82  os << defaultMean << "\n";
83  os.precision(pr);
84  return os;
85 #endif
86 }
87 
88 std::istream & RandExpZiggurat::get ( std::istream & is ) {
89  std::string inName;
90  is >> inName;
91  if (inName != name()) {
92  is.clear(std::ios::badbit | is.rdstate());
93  std::cerr << "Mismatch when expecting to read state of a "
94  << name() << " distribution\n"
95  << "Name found was " << inName
96  << "\nistream is left in the badbit state\n";
97  return is;
98  }
99  if (possibleKeywordInput(is, "Uvec", defaultMean)) {
100  std::vector<unsigned long> t(2);
101  is >> defaultMean >> t[0] >> t[1]; defaultMean = DoubConv::longs2double(t);
102  return is;
103  }
104  // is >> defaultMean encompassed by possibleKeywordInput
105  return is;
106 }
107 
108 
109 float RandExpZiggurat::ziggurat_efix(unsigned long jz,HepRandomEngine* anEngine)
110 {
112 
113  unsigned long iz=jz&255;
114 
115  float x;
116  for(;;)
117  {
118  if(iz==0) return (7.69711-std::log(ziggurat_UNI(anEngine))); /* iz==0 */
119  x=jz*we[iz];
120  if( fe[iz]+ziggurat_UNI(anEngine)*(fe[iz-1]-fe[iz]) < std::exp(-x) ) return (x);
121 
122  /* initiate, try to exit for(;;) loop */
123  jz=ziggurat_SHR3(anEngine);
124  iz=(jz&255);
125  if(jz<ke[iz]) return (jz*we[iz]);
126  }
127 }
128 
130 {
131  const double rzm1 = 2147483648.0, rzm2 = 4294967296.;
132  double dn=3.442619855899,tn=dn,vn=9.91256303526217e-3, q;
133  double de=7.697117470131487, te=de, ve=3.949659822581572e-3;
134  int i;
135 
136 /* Set up tables for RNOR */
137  q=vn/std::exp(-.5*dn*dn);
138  kn[0]=(unsigned long)((dn/q)*rzm1);
139  kn[1]=0;
140 
141  wn[0]=q/rzm1;
142  wn[127]=dn/rzm1;
143 
144  fn[0]=1.;
145  fn[127]=std::exp(-.5*dn*dn);
146 
147  for(i=126;i>=1;i--) {
148  dn=std::sqrt(-2.*std::log(vn/dn+std::exp(-.5*dn*dn)));
149  kn[i+1]=(unsigned long)((dn/tn)*rzm1);
150  tn=dn;
151  fn[i]=std::exp(-.5*dn*dn);
152  wn[i]=dn/rzm1;
153  }
154 
155 /* Set up tables for REXP */
156  q = ve/std::exp(-de);
157  ke[0]=(unsigned long)((de/q)*rzm2);
158  ke[1]=0;
159 
160  we[0]=q/rzm2;
161  we[255]=de/rzm2;
162 
163  fe[0]=1.;
164  fe[255]=std::exp(-de);
165 
166  for(i=254;i>=1;i--) {
167  de=-std::log(ve/de+std::exp(-de));
168  ke[i+1]= (unsigned long)((de/te)*rzm2);
169  te=de;
170  fe[i]=std::exp(-de);
171  we[i]=de/rzm2;
172  }
173  ziggurat_is_init=true;
174  return true;
175 }
176 
177 } // namespace CLHEP