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