Geant4  10.00.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 bool RandExpZiggurat::ziggurat_is_init=RandExpZiggurat::ziggurat_init();
11 unsigned long RandExpZiggurat::kn[128], RandExpZiggurat::ke[256];
12 float RandExpZiggurat::wn[128],RandExpZiggurat::fn[128],RandExpZiggurat::we[256],RandExpZiggurat::fe[256];
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  if(!ziggurat_is_init) ziggurat_init();
25 }
26 
27 double RandExpZiggurat::operator()()
28 {
29  return fire( defaultMean );
30 }
31 
32 void RandExpZiggurat::shootArray( const int size, float* vect, float mean )
33 {
34  for (int i=0; i<size; ++i) vect[i] = shoot(mean);
35 }
36 
37 void RandExpZiggurat::shootArray( const int size, double* vect, double mean )
38 {
39  for (int i=0; i<size; ++i) vect[i] = shoot(mean);
40 }
41 
42 void RandExpZiggurat::shootArray(HepRandomEngine* anEngine, const int size, float* vect, float mean )
43 {
44  for (int i=0; i<size; ++i) vect[i] = shoot(anEngine, mean);
45 }
46 
47 void RandExpZiggurat::shootArray(HepRandomEngine* anEngine, const int size, double* vect, double mean )
48 {
49  for (int i=0; i<size; ++i) vect[i] = shoot(anEngine, mean);
50 }
51 
52 void RandExpZiggurat::fireArray( const int size, float* vect)
53 {
54  for (int i=0; i<size; ++i) vect[i] = fire( defaultMean );
55 }
56 
57 void RandExpZiggurat::fireArray( const int size, double* vect)
58 {
59  for (int i=0; i<size; ++i) vect[i] = fire( defaultMean );
60 }
61 
62 void RandExpZiggurat::fireArray( const int size, float* vect, float mean )
63 {
64  for (int i=0; i<size; ++i) vect[i] = fire( mean );
65 }
66 
67 void RandExpZiggurat::fireArray( const int size, double* vect, double mean )
68 {
69  for (int i=0; i<size; ++i) vect[i] = fire( mean );
70 }
71 
72 std::ostream & RandExpZiggurat::put ( std::ostream & os ) const {
73  int pr=os.precision(20);
74  std::vector<unsigned long> t(2);
75  os << " " << name() << "\n";
76  os << "Uvec" << "\n";
77  t = DoubConv::dto2longs(defaultMean);
78  os << defaultMean << " " << t[0] << " " << t[1] << "\n";
79  os.precision(pr);
80  return os;
81 #ifdef REMOVED
82  int pr=os.precision(20);
83  os << " " << name() << "\n";
84  os << defaultMean << "\n";
85  os.precision(pr);
86  return os;
87 #endif
88 }
89 
90 std::istream & RandExpZiggurat::get ( std::istream & is ) {
91  std::string inName;
92  is >> inName;
93  if (inName != name()) {
94  is.clear(std::ios::badbit | is.rdstate());
95  std::cerr << "Mismatch when expecting to read state of a "
96  << name() << " distribution\n"
97  << "Name found was " << inName
98  << "\nistream is left in the badbit state\n";
99  return is;
100  }
101  if (possibleKeywordInput(is, "Uvec", defaultMean)) {
102  std::vector<unsigned long> t(2);
103  is >> defaultMean >> t[0] >> t[1]; defaultMean = DoubConv::longs2double(t);
104  return is;
105  }
106  // is >> defaultMean encompassed by possibleKeywordInput
107  return is;
108 }
109 
110 
111 float RandExpZiggurat::ziggurat_efix(unsigned long jz,HepRandomEngine* anEngine)
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 
129 bool RandExpZiggurat::ziggurat_init()
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
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