Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
RandGauss.cc
Go to the documentation of this file.
1// -*- C++ -*-
2//
3// -----------------------------------------------------------------------
4// HEP Random
5// --- RandGauss ---
6// class implementation file
7// -----------------------------------------------------------------------
8// This file is part of Geant4 (simulation toolkit for HEP).
9
10// =======================================================================
11// Gabriele Cosmo - Created: 5th September 1995
12// - Added methods to shoot arrays: 28th July 1997
13// J.Marraffino - Added default arguments as attributes and
14// operator() with arguments. Introduced method normal()
15// for computation in fire(): 16th Feb 1998
16// Gabriele Cosmo - Relocated static data from HepRandom: 5th Jan 1999
17// M Fischler - Copy constructor should supply right engine to HepRandom:
18// 1/26/00.
19// M Fischler - Workaround for problem of non-reproducing saveEngineStatus
20// by saving cached gaussian. March 2000.
21// M Fischler - Avoiding hang when file not found in restoreEngineStatus
22// 12/3/04
23// M Fischler - put and get to/from streams 12/8/04
24// M Fischler - save and restore dist to streams 12/20/04
25// M Fischler - put/get to/from streams uses pairs of ulongs when
26// storing doubles avoid problems with precision.
27// Similarly for saveEngineStatus and RestoreEngineStatus
28// and for save/restore distState
29// Care was taken that old-form output can still be read back.
30// 4/14/05
31//
32// =======================================================================
33
36#include <cmath> // for std::log()
37#include <iostream>
38#include <string.h> // for strcmp
39#include <string>
40#include <vector>
41
42namespace CLHEP {
43
44std::string RandGauss::name() const {return "RandGauss";}
46
47// Initialisation of static data
48CLHEP_THREAD_LOCAL bool RandGauss::set_st = false;
49CLHEP_THREAD_LOCAL double RandGauss::nextGauss_st = 0.0;
50
53
57
58double RandGauss::operator()( double mean, double stdDev ) {
59 return fire( mean, stdDev );
60}
61
63{
64 // Gaussian random numbers are generated two at the time, so every other
65 // time this is called we just return a number generated the time before.
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;
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
91void RandGauss::shootArray( const int size, double* vect,
92 double mean, double stdDev )
93{
94 for( double* v = vect; v != vect + size; ++v )
95 *v = shoot(mean,stdDev);
96}
97
99{
100 // Gaussian random numbers are generated two at the time, so every other
101 // time this is called we just return a number generated the time before.
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
125 const int size, double* vect,
126 double mean, double stdDev )
127{
128 for( double* v = vect; v != vect + size; ++v )
129 *v = shoot(anEngine,mean,stdDev);
130}
131
133{
134 // Gaussian random numbers are generated two at the time, so every other
135 // time this is called we just return a number generated the time before.
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
158void RandGauss::fireArray( const int size, double* vect)
159{
160 for( double* v = vect; v != vect + size; ++v )
162}
163
164void RandGauss::fireArray( const int size, double* vect,
165 double mean, double stdDev )
166{
167 for( double* v = vect; v != vect + size; ++v )
168 *v = fire( mean, stdDev );
169}
170
172{
173 return set_st;
174}
175
176void RandGauss::setFlag( bool val )
177{
178 set_st = val;
179}
180
182{
183 return nextGauss_st;
184}
185
186void RandGauss::setVal( double nextVal )
187{
188 nextGauss_st = nextVal;
189}
190
191void RandGauss::saveEngineStatus ( const char filename[] ) {
192
193 // First save the engine status just like the base class would do:
194 getTheEngine()->saveStatus( filename );
195
196 // Now append the cached variate, if any:
197
198 std::ofstream outfile ( filename, std::ios::app );
199
200 if ( getFlag() ) {
201 std::vector<unsigned long> t(2);
203 outfile << "RANDGAUSS CACHED_GAUSSIAN: Uvec "
204 << getVal() << " " << t[0] << " " << t[1] << "\n";
205 } else {
206 outfile << "RANDGAUSS NO_CACHED_GAUSSIAN: 0 \n" ;
207 }
208
209} // saveEngineStatus
210
211void RandGauss::restoreEngineStatus( const char filename[] ) {
212
213 // First restore the engine status just like the base class would do:
214 getTheEngine()->restoreStatus( filename );
215
216 // Now find the line describing the cached variate:
217
218 std::ifstream infile ( filename, std::ios::in );
219 if (!infile) return;
220
221 char inputword[] = "NO_KEYWORD "; // leaves room for 14 characters plus \0
222 while (true) {
223 infile.width(13);
224 infile >> inputword;
225 if (strcmp(inputword,"RANDGAUSS")==0) break;
226 if (infile.eof()) break;
227 // If the file ends without the RANDGAUSS line, that means this
228 // was a file produced by an earlier version of RandGauss. We will
229 // replicated the old behavior in that case: set_st is cleared.
230 }
231
232 // Then read and use the caching info:
233
234 if (strcmp(inputword,"RANDGAUSS")==0) {
235 char setword[40]; // the longest, staticFirstUnusedBit: has length 21
236 infile.width(39);
237 infile >> setword; // setword should be CACHED_GAUSSIAN:
238 if (strcmp(setword,"CACHED_GAUSSIAN:") ==0) {
239 if (possibleKeywordInput(infile, "Uvec", nextGauss_st)) {
240 std::vector<unsigned long> t(2);
241 infile >> nextGauss_st >> t[0] >> t[1];
242 nextGauss_st = DoubConv::longs2double(t);
243 }
244 // is >> nextGauss_st encompassed by possibleKeywordInput
245 setFlag(true);
246 } else {
247 setFlag(false);
248 infile >> nextGauss_st; // because a 0 will have been output
249 }
250 } else {
251 setFlag(false);
252 }
253
254} // restoreEngineStatus
255
256 // Save and restore to/from streams
257
258std::ostream & RandGauss::put ( std::ostream & os ) const {
259 os << name() << "\n";
260 long prec = os.precision(20);
261 std::vector<unsigned long> t(2);
262 os << "Uvec\n";
264 os << defaultMean << " " << t[0] << " " << t[1] << "\n";
266 os << defaultStdDev << " " << t[0] << " " << t[1] << "\n";
267 if ( set ) {
268 t = DoubConv::dto2longs(nextGauss);
269 os << "nextGauss " << nextGauss << " " << t[0] << " " << t[1] << "\n";
270 } else {
271 os << "no_cached_nextGauss \n";
272 }
273 os.precision(prec);
274 return os;
275} // put
276
277std::istream & RandGauss::get ( std::istream & is ) {
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 read state of a "
283 << name() << " distribution\n"
284 << "Name found was " << inName
285 << "\nistream is left in the badbit state\n";
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]; defaultMean = DoubConv::longs2double(t);
293 is >> defaultStdDev>>t[0]>>t[1]; defaultStdDev = DoubConv::longs2double(t);
294 std::string ng;
295 is >> ng;
296 set = false;
297 if (ng == "nextGauss") {
298 is >> nextGauss >> t[0] >> t[1]; nextGauss = DoubConv::longs2double(t);
299 set = true;
300 }
301 return is;
302 }
303 // is >> c1 encompassed by possibleKeywordInput
304 is >> defaultMean >> c2 >> defaultStdDev;
305 if ( (!is) || (c1 != "Mean:") || (c2 != "Sigma:") ) {
306 std::cerr << "i/o problem while expecting to read state of a "
307 << name() << " distribution\n"
308 << "default mean and/or sigma could not be read\n";
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 state of RandGauss\n";
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 keyword of RandGauss:" << c2
324 << "\nistream is left in the badbit state\n";
325 }
326 return is;
327} // get
328
329 // Static save and restore to/from streams
330
331std::ostream & RandGauss::saveDistState ( std::ostream & os ) {
332 long prec = os.precision(20);
333 std::vector<unsigned long> t(2);
334 os << distributionName() << "\n";
335 os << "Uvec\n";
336 if ( getFlag() ) {
338 os << "nextGauss_st " << getVal() << " " << t[0] << " " << t[1] << "\n";
339 } else {
340 os << "no_cached_nextGauss_st \n";
341 }
342 os.precision(prec);
343 return os;
344}
345
346std::istream & RandGauss::restoreDistState ( std::istream & is ) {
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 read static state of a "
352 << distributionName() << " distribution\n"
353 << "Name found was " << inName
354 << "\nistream is left in the badbit state\n";
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 possibleKeywordInput
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 state of static RandGauss\n";
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 keyword of static RandGauss:" << c2
385 << "\nistream is left in the badbit state\n";
386 }
387 return is;
388}
389
390std::ostream & RandGauss::saveFullState ( std::ostream & os ) {
392 saveDistState(os);
393 return os;
394}
395
396std::istream & RandGauss::restoreFullState ( std::istream & is ) {
399 return is;
400}
401
402} // namespace CLHEP
403
static double longs2double(const std::vector< unsigned long > &v)
Definition DoubConv.cc:110
static std::vector< unsigned long > dto2longs(double d)
Definition DoubConv.cc:94
virtual void restoreStatus(const char filename[]="Config.conf")=0
virtual double flat()=0
virtual void saveStatus(const char filename[]="Config.conf") const =0
static HepRandomEngine * getTheEngine()
Definition Random.cc:268
static std::ostream & saveFullState(std::ostream &os)
Definition Random.cc:288
static std::istream & restoreFullState(std::istream &is)
Definition Random.cc:293
static std::ostream & saveDistState(std::ostream &os)
Definition RandGauss.cc:331
std::string name() const
Definition RandGauss.cc:44
static std::istream & restoreFullState(std::istream &is)
Definition RandGauss.cc:396
static void restoreEngineStatus(const char filename[]="Config.conf")
Definition RandGauss.cc:211
static double shoot()
Definition RandGauss.cc:62
static bool getFlag()
Definition RandGauss.cc:171
static double getVal()
Definition RandGauss.cc:181
std::istream & get(std::istream &is)
Definition RandGauss.cc:277
double defaultStdDev
Definition RandGauss.h:152
std::ostream & put(std::ostream &os) const
Definition RandGauss.cc:258
void fireArray(const int size, double *vect)
Definition RandGauss.cc:158
static void shootArray(const int size, double *vect, double mean=0.0, double stdDev=1.0)
Definition RandGauss.cc:91
static void setVal(double nextVal)
Definition RandGauss.cc:186
HepRandomEngine & engine()
Definition RandGauss.cc:45
static std::ostream & saveFullState(std::ostream &os)
Definition RandGauss.cc:390
static void setFlag(bool val)
Definition RandGauss.cc:176
static std::string distributionName()
Definition RandGauss.h:99
std::shared_ptr< HepRandomEngine > localEngine
Definition RandGauss.h:154
virtual ~RandGauss()
Definition RandGauss.cc:51
static std::istream & restoreDistState(std::istream &is)
Definition RandGauss.cc:346
virtual double operator()()
Definition RandGauss.cc:54
static void saveEngineStatus(const char filename[]="Config.conf")
Definition RandGauss.cc:191
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
#define CLHEP_THREAD_LOCAL