CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

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