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 )
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);
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 ) {
391  return is;
392 }
393 
394 } // namespace CLHEP
395 
CLHEP::HepRandom::saveFullState
static std::ostream & saveFullState(std::ostream &os)
Definition: Random.cc:186
CLHEP::RandGauss::engine
HepRandomEngine & engine()
Definition: RandGauss.cc:44
CLHEP::HepRandomEngine
Definition: Matrix/CLHEP/Random/RandomEngine.h:55
CLHEP::RandGauss::restoreDistState
static std::istream & restoreDistState(std::istream &is)
Definition: RandGauss.cc:338
is
HepRotation and so forth isNear() norm2() rectify() static Rotation row1 row4(To avoid bloat in the code pulled in for programs which don 't use all these features, we split the implementation .cc files. Only isNear() goes into the original Rotation.cc) --------------------------------------- HepAxisAngle and HepEulerAngles classes --------------------------------------- These classes are very useful and simple structures for holding the result of a nice intuituve decomposition of a rotation there is no longer much content in the distinct ZOOM PhysicsVectors library The only content left in the library is the object files representing the various Exception objects When we build the CLHEP classes for the ZOOM we will set up so as to use ZOOM SpaceVector is(but we can disable namespace usage and most of our users do so at this point). What I do is leave Hep3Vector in the global namespace
CLHEP::HepRandom::restoreFullState
static std::istream & restoreFullState(std::istream &is)
Definition: Random.cc:191
CLHEP::RandGauss::restoreFullState
static std::istream & restoreFullState(std::istream &is)
Definition: RandGauss.cc:388
CLHEP::RandGauss::normal
double normal()
Definition: RandGauss.cc:131
CLHEP::RandGauss::getVal
static double getVal()
Definition: Matrix/CLHEP/Random/RandGauss.h:146
CLHEP::DoubConv::longs2double
static double longs2double(const std::vector< unsigned long > &v)
Definition: DoubConv.cc:106
CLHEP::RandGauss::operator()
virtual double operator()()
Definition: RandGauss.cc:53
CLHEP::RandGauss::getFlag
static bool getFlag()
Definition: Matrix/CLHEP/Random/RandGauss.h:112
CLHEP::RandGauss::setFlag
static void setFlag(bool val)
Definition: Matrix/CLHEP/Random/RandGauss.h:114
CLHEP::RandGauss::get
std::istream & get(std::istream &is)
Definition: RandGauss.cc:262
CLHEP::RandGauss::~RandGauss
virtual ~RandGauss()
Definition: RandGauss.cc:50
size
user code seldom needs to call this function directly ZMerrno whether or not they are still recorded ZMerrno size() Return the(integer) number of ZMthrow 'n exceptions currently recorded. 5) ZMerrno.clear() Set an internal counter to zero. This counter is available(see next function) to user code to track ZMthrow 'n exceptions that have occurred during any arbitrary time interval. 6) ZMerrno.countSinceCleared() Return the(integer) number of ZMthrow 'n exceptions that have been recorded via ZMerrno.write()
CLHEP::RandGauss::shoot
static double shoot()
Definition: RandGauss.cc:61
CLHEP::RandGauss::defaultStdDev
double defaultStdDev
Definition: Matrix/CLHEP/Random/RandGauss.h:153
CLHEP::RandGauss::saveDistState
static std::ostream & saveDistState(std::ostream &os)
Definition: RandGauss.cc:323
CLHEP::RandGauss::fireArray
void fireArray(const int size, double *vect)
Definition: RandGauss.cc:157
CLHEP
Definition: ClhepVersion.h:13
CLHEP::RandGauss::restoreEngineStatus
static void restoreEngineStatus(const char filename[]="Config.conf")
Definition: RandGauss.cc:190
v
they are gone ZOOM Features Discontinued The following features of the ZOOM package were felt to be extreme overkill These have been after checking that no existing user code was utilizing as in SpaceVector v
Definition: keyMergeIssues.doc:324
CLHEP::HepRandomEngine::saveStatus
virtual void saveStatus(const char filename[]="Config.conf") const =0
CLHEP::RandGauss::defaultMean
double defaultMean
Definition: Matrix/CLHEP/Random/RandGauss.h:152
CLHEP::RandGauss::distributionName
static std::string distributionName()
Definition: Matrix/CLHEP/Random/RandGauss.h:100
CLHEP::possibleKeywordInput
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
Definition: Matrix/CLHEP/Random/RandomEngine.h:168
CLHEP::HepRandomEngine::restoreStatus
virtual void restoreStatus(const char filename[]="Config.conf")=0
CLHEP::RandGauss::name
std::string name() const
Definition: RandGauss.cc:43
CLHEP::HepRandomEngine::flat
virtual double flat()=0
CLHEP::RandGauss::saveEngineStatus
static void saveEngineStatus(const char filename[]="Config.conf")
Definition: RandGauss.cc:170
CLHEP::RandGauss::put
std::ostream & put(std::ostream &os) const
Definition: RandGauss.cc:237
CLHEP::DoubConv::dto2longs
static std::vector< unsigned long > dto2longs(double d)
Definition: DoubConv.cc:90
CLHEP::RandGauss::saveFullState
static std::ostream & saveFullState(std::ostream &os)
Definition: RandGauss.cc:382
CLHEP::RandGauss::localEngine
shared_ptr< HepRandomEngine > localEngine
Definition: Matrix/CLHEP/Random/RandGauss.h:155
in
it has advantages For I leave the ZMthrows in
Definition: keyMergeIssues.doc:62
CLHEP::HepRandom::getTheEngine
static HepRandomEngine * getTheEngine()
Definition: Random.cc:166
x
any side effects of that construction would occur twice The semantics of throw x
Definition: whyZMthrowRethrows.txt:37
CLHEP::RandGauss::shootArray
static void shootArray(const int size, double *vect, double mean=0.0, double stdDev=1.0)
Definition: RandGauss.cc:90
CLHEP::RandGauss::setVal
static void setVal(double nextVal)
Definition: Matrix/CLHEP/Random/RandGauss.h:148
CLHEP::RandGauss::fire
double fire()