| Geant4 Cross Reference |
1 // -*- C++ -*- 1
2 //
3 // -------------------------------------------
4 // HEP Random
5 // --- RandGauss ---
6 // class implementation f
7 // -------------------------------------------
8 // This file is part of Geant4 (simulation too
9
10 // ===========================================
11 // Gabriele Cosmo - Created: 5th September 199
12 // - Added methods to shoot arr
13 // J.Marraffino - Added default arguments as
14 // operator() with arguments.
15 // for computation in fire():
16 // Gabriele Cosmo - Relocated static data from
17 // M Fischler - Copy constructor should su
18 // 1/26/00.
19 // M Fischler - Workaround for problem of
20 // by saving cached gaussian. March 20
21 // M Fischler - Avoiding hang when file no
22 // 12/3/04
23 // M Fischler - put and get to/from stream
24 // M Fischler - save and restore dist to s
25 // M Fischler - put/get to/from streams uses
26 // storing doubles avoid problems with
27 // Similarly for saveEngineStatus and R
28 // and for save/restore distState
29 // Care was taken that old-form output
30 // 4/14/05
31 //
32 // ===========================================
33
34 #include "CLHEP/Random/RandGauss.h"
35 #include "CLHEP/Random/DoubConv.h"
36 #include <cmath> // for std::log()
37 #include <iostream>
38 #include <string.h> // for strcmp
39 #include <string>
40 #include <vector>
41
42 namespace CLHEP {
43
44 std::string RandGauss::name() const {return "R
45 HepRandomEngine & RandGauss::engine() {return
46
47 // Initialisation of static data
48 CLHEP_THREAD_LOCAL bool RandGauss::set_st = fa
49 CLHEP_THREAD_LOCAL double RandGauss::nextGauss
50
51 RandGauss::~RandGauss() {
52 }
53
54 double RandGauss::operator()() {
55 return fire( defaultMean, defaultStdDev );
56 }
57
58 double RandGauss::operator()( double mean, dou
59 return fire( mean, stdDev );
60 }
61
62 double RandGauss::shoot()
63 {
64 // Gaussian random numbers are generated two
65 // time this is called we just return a numb
66
67 if ( getFlag() ) {
68 setFlag(false);
69 double x = getVal();
70 return x;
71 // return getVal();
72 }
73
74 double r;
75 double v1,v2,fac,val;
76 HepRandomEngine* anEngine = HepRandom::getTh
77
78 do {
79 v1 = 2.0 * anEngine->flat() - 1.0;
80 v2 = 2.0 * anEngine->flat() - 1.0;
81 r = v1*v1 + v2*v2;
82 } while ( r > 1.0 );
83
84 fac = std::sqrt(-2.0*std::log(r)/r);
85 val = v1*fac;
86 setVal(val);
87 setFlag(true);
88 return v2*fac;
89 }
90
91 void RandGauss::shootArray( const int size, do
92 double mean, doubl
93 {
94 for( double* v = vect; v != vect + size; ++v
95 *v = shoot(mean,stdDev);
96 }
97
98 double RandGauss::shoot( HepRandomEngine* anEn
99 {
100 // Gaussian random numbers are generated two
101 // time this is called we just return a numb
102
103 if ( getFlag() ) {
104 setFlag(false);
105 return getVal();
106 }
107
108 double r;
109 double v1,v2,fac,val;
110
111 do {
112 v1 = 2.0 * anEngine->flat() - 1.0;
113 v2 = 2.0 * anEngine->flat() - 1.0;
114 r = v1*v1 + v2*v2;
115 } while ( r > 1.0 );
116
117 fac = std::sqrt( -2.0*std::log(r)/r);
118 val = v1*fac;
119 setVal(val);
120 setFlag(true);
121 return v2*fac;
122 }
123
124 void RandGauss::shootArray( HepRandomEngine* a
125 const int size, do
126 double mean, doubl
127 {
128 for( double* v = vect; v != vect + size; ++v
129 *v = shoot(anEngine,mean,stdDev);
130 }
131
132 double RandGauss::normal()
133 {
134 // Gaussian random numbers are generated two
135 // time this is called we just return a numb
136
137 if ( set ) {
138 set = false;
139 return nextGauss;
140 }
141
142 double r;
143 double v1,v2,fac,val;
144
145 do {
146 v1 = 2.0 * localEngine->flat() - 1.0;
147 v2 = 2.0 * localEngine->flat() - 1.0;
148 r = v1*v1 + v2*v2;
149 } while ( r > 1.0 );
150
151 fac = std::sqrt(-2.0*std::log(r)/r);
152 val = v1*fac;
153 nextGauss = val;
154 set = true;
155 return v2*fac;
156 }
157
158 void RandGauss::fireArray( const int size, dou
159 {
160 for( double* v = vect; v != vect + size; ++v
161 *v = fire( defaultMean, defaultStdDev );
162 }
163
164 void RandGauss::fireArray( const int size, dou
165 double mean, double
166 {
167 for( double* v = vect; v != vect + size; ++v
168 *v = fire( mean, stdDev );
169 }
170
171 bool RandGauss::getFlag()
172 {
173 return set_st;
174 }
175
176 void RandGauss::setFlag( bool val )
177 {
178 set_st = val;
179 }
180
181 double RandGauss::getVal()
182 {
183 return nextGauss_st;
184 }
185
186 void RandGauss::setVal( double nextVal )
187 {
188 nextGauss_st = nextVal;
189 }
190
191 void RandGauss::saveEngineStatus ( const char
192
193 // First save the engine status just like th
194 getTheEngine()->saveStatus( filename );
195
196 // Now append the cached variate, if any:
197
198 std::ofstream outfile ( filename, std::ios::
199
200 if ( getFlag() ) {
201 std::vector<unsigned long> t(2);
202 t = DoubConv::dto2longs(getVal());
203 outfile << "RANDGAUSS CACHED_GAUSSIAN: Uve
204 << getVal() << " " << t[0] << " " << t[1
205 } else {
206 outfile << "RANDGAUSS NO_CACHED_GAUSSIAN:
207 }
208
209 } // saveEngineStatus
210
211 void RandGauss::restoreEngineStatus( const cha
212
213 // First restore the engine status just like
214 getTheEngine()->restoreStatus( filename );
215
216 // Now find the line describing the cached v
217
218 std::ifstream infile ( filename, std::ios::i
219 if (!infile) return;
220
221 char inputword[] = "NO_KEYWORD "; // leav
222 while (true) {
223 infile.width(13);
224 infile >> inputword;
225 if (strcmp(inputword,"RANDGAUSS")==0) brea
226 if (infile.eof()) break;
227 // If the file ends without the RANDGAUSS li
228 // was a file produced by an earlier version
229 // replicated the old behavior in that case:
230 }
231
232 // Then read and use the caching info:
233
234 if (strcmp(inputword,"RANDGAUSS")==0) {
235 char setword[40]; // the longest, staticFi
236 infile.width(39);
237 infile >> setword; // setword should be C
238 if (strcmp(setword,"CACHED_GAUSSIAN:") ==0
239 if (possibleKeywordInput(infile, "Uvec",
240 std::vector<unsigned long> t(2);
241 infile >> nextGauss_st >> t[0] >> t[1]
242 nextGauss_st = DoubConv::longs2double(
243 }
244 // is >> nextGauss_st encompassed by pos
245 setFlag(true);
246 } else {
247 setFlag(false);
248 infile >> nextGauss_st; // because a 0 w
249 }
250 } else {
251 setFlag(false);
252 }
253
254 } // restoreEngineStatus
255
256 // Save and restore to/from streams
257
258 std::ostream & RandGauss::put ( std::ostream &
259 os << name() << "\n";
260 long prec = os.precision(20);
261 std::vector<unsigned long> t(2);
262 os << "Uvec\n";
263 t = DoubConv::dto2longs(defaultMean);
264 os << defaultMean << " " << t[0] << " " << t
265 t = DoubConv::dto2longs(defaultStdDev);
266 os << defaultStdDev << " " << t[0] << " " <<
267 if ( set ) {
268 t = DoubConv::dto2longs(nextGauss);
269 os << "nextGauss " << nextGauss << " " <<
270 } else {
271 os << "no_cached_nextGauss \n";
272 }
273 os.precision(prec);
274 return os;
275 } // put
276
277 std::istream & RandGauss::get ( std::istream &
278 std::string inName;
279 is >> inName;
280 if (inName != name()) {
281 is.clear(std::ios::badbit | is.rdstate());
282 std::cerr << "Mismatch when expecting to r
283 << name() << " distribution\n"
284 << "Name found was " << inName
285 << "\nistream is left in the badbit st
286 return is;
287 }
288 std::string c1;
289 std::string c2;
290 if (possibleKeywordInput(is, "Uvec", c1)) {
291 std::vector<unsigned long> t(2);
292 is >> defaultMean >> t[0] >> t[1]; default
293 is >> defaultStdDev>>t[0]>>t[1]; defaultSt
294 std::string ng;
295 is >> ng;
296 set = false;
297 if (ng == "nextGauss") {
298 is >> nextGauss >> t[0] >> t[1]; nextGau
299 set = true;
300 }
301 return is;
302 }
303 // is >> c1 encompassed by possibleKeywordIn
304 is >> defaultMean >> c2 >> defaultStdDev;
305 if ( (!is) || (c1 != "Mean:") || (c2 != "Sig
306 std::cerr << "i/o problem while expecting
307 << name() << " distribution\n"
308 << "default mean and/or sigma could no
309 return is;
310 }
311 is >> c1 >> c2 >> nextGauss;
312 if ( (!is) || (c1 != "RANDGAUSS") ) {
313 is.clear(std::ios::badbit | is.rdstate());
314 std::cerr << "Failure when reading caching
315 return is;
316 }
317 if (c2 == "CACHED_GAUSSIAN:") {
318 set = true;
319 } else if (c2 == "NO_CACHED_GAUSSIAN:") {
320 set = false;
321 } else {
322 is.clear(std::ios::badbit | is.rdstate());
323 std::cerr << "Unexpected caching state key
324 << "\nistream is left in the badbit st
325 }
326 return is;
327 } // get
328
329 // Static save and restore to/from streams
330
331 std::ostream & RandGauss::saveDistState ( std:
332 long prec = os.precision(20);
333 std::vector<unsigned long> t(2);
334 os << distributionName() << "\n";
335 os << "Uvec\n";
336 if ( getFlag() ) {
337 t = DoubConv::dto2longs(getVal());
338 os << "nextGauss_st " << getVal() << " " <
339 } else {
340 os << "no_cached_nextGauss_st \n";
341 }
342 os.precision(prec);
343 return os;
344 }
345
346 std::istream & RandGauss::restoreDistState ( s
347 std::string inName;
348 is >> inName;
349 if (inName != distributionName()) {
350 is.clear(std::ios::badbit | is.rdstate());
351 std::cerr << "Mismatch when expecting to r
352 << distributionName() << " distrib
353 << "Name found was " << inName
354 << "\nistream is left in the badbit st
355 return is;
356 }
357 std::string c1;
358 std::string c2;
359 if (possibleKeywordInput(is, "Uvec", c1)) {
360 std::vector<unsigned long> t(2);
361 std::string ng;
362 is >> ng;
363 setFlag (false);
364 if (ng == "nextGauss_st") {
365 is >> nextGauss_st >> t[0] >> t[1];
366 nextGauss_st = DoubConv::longs2double(t)
367 setFlag (true);
368 }
369 return is;
370 }
371 // is >> c1 encompassed by possibleKeywordIn
372 is >> c2 >> nextGauss_st;
373 if ( (!is) || (c1 != "RANDGAUSS") ) {
374 is.clear(std::ios::badbit | is.rdstate());
375 std::cerr << "Failure when reading caching
376 return is;
377 }
378 if (c2 == "CACHED_GAUSSIAN:") {
379 setFlag(true);
380 } else if (c2 == "NO_CACHED_GAUSSIAN:") {
381 setFlag(false);
382 } else {
383 is.clear(std::ios::badbit | is.rdstate());
384 std::cerr << "Unexpected caching state key
385 << "\nistream is left in the badbit st
386 }
387 return is;
388 }
389
390 std::ostream & RandGauss::saveFullState ( std:
391 HepRandom::saveFullState(os);
392 saveDistState(os);
393 return os;
394 }
395
396 std::istream & RandGauss::restoreFullState ( s
397 HepRandom::restoreFullState(is);
398 restoreDistState(is);
399 return is;
400 }
401
402 } // namespace CLHEP
403
404