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 
41 namespace CLHEP {
42 
43 std::string RandGauss::name() const {return "RandGauss";}
45 
46 // Initialisation of static data
47 bool RandGauss::set_st = false;
48 double RandGauss::nextGauss_st = 0.0;
49 
51 }
52 
54  return fire( defaultMean, defaultStdDev );
55 }
56 
57 double 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 
90 void 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 
97 double RandGauss::shoot( HepRandomEngine* anEngine )
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 
157 void RandGauss::fireArray( const int size, double* vect)
158 {
159  for( double* v = vect; v != vect + size; ++v )
160  *v = fire( defaultMean, defaultStdDev );
161 }
162 
163 void 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 
170 void 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 
190 void 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 
237 std::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 
262 std::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 
323 std::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 
338 std::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 
382 std::ostream & RandGauss::saveFullState ( std::ostream & os ) {
384  saveDistState(os);
385  return os;
386 }
387 
388 std::istream & RandGauss::restoreFullState ( std::istream & is ) {
390  restoreDistState(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)