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