Geant4  10.00.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 bool RandGauss::set_st = false;
47 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 void RandGauss::saveEngineStatus ( const char filename[] ) {
170 
171  // First save the engine status just like the base class would do:
172  getTheEngine()->saveStatus( filename );
173 
174  // Now append the cached variate, if any:
175 
176  std::ofstream outfile ( filename, std::ios::app );
177 
178  if ( getFlag() ) {
179  std::vector<unsigned long> t(2);
180  t = DoubConv::dto2longs(getVal());
181  outfile << "RANDGAUSS CACHED_GAUSSIAN: Uvec "
182  << getVal() << " " << t[0] << " " << t[1] << "\n";
183  } else {
184  outfile << "RANDGAUSS NO_CACHED_GAUSSIAN: 0 \n" ;
185  }
186 
187 } // saveEngineStatus
188 
189 void RandGauss::restoreEngineStatus( const char filename[] ) {
190 
191  // First restore the engine status just like the base class would do:
192  getTheEngine()->restoreStatus( filename );
193 
194  // Now find the line describing the cached variate:
195 
196  std::ifstream infile ( filename, std::ios::in );
197  if (!infile) return;
198 
199  char inputword[] = "NO_KEYWORD "; // leaves room for 14 characters plus \0
200  while (true) {
201  infile.width(13);
202  infile >> inputword;
203  if (strcmp(inputword,"RANDGAUSS")==0) break;
204  if (infile.eof()) break;
205  // If the file ends without the RANDGAUSS line, that means this
206  // was a file produced by an earlier version of RandGauss. We will
207  // replicated the old behavior in that case: set_st is cleared.
208  }
209 
210  // Then read and use the caching info:
211 
212  if (strcmp(inputword,"RANDGAUSS")==0) {
213  char setword[40]; // the longest, staticFirstUnusedBit: has length 21
214  infile.width(39);
215  infile >> setword; // setword should be CACHED_GAUSSIAN:
216  if (strcmp(setword,"CACHED_GAUSSIAN:") ==0) {
217  if (possibleKeywordInput(infile, "Uvec", nextGauss_st)) {
218  std::vector<unsigned long> t(2);
219  infile >> nextGauss_st >> t[0] >> t[1];
220  nextGauss_st = DoubConv::longs2double(t);
221  }
222  // is >> nextGauss_st encompassed by possibleKeywordInput
223  setFlag(true);
224  } else {
225  setFlag(false);
226  infile >> nextGauss_st; // because a 0 will have been output
227  }
228  } else {
229  setFlag(false);
230  }
231 
232 } // restoreEngineStatus
233 
234  // Save and restore to/from streams
235 
236 std::ostream & RandGauss::put ( std::ostream & os ) const {
237  os << name() << "\n";
238  int prec = os.precision(20);
239  std::vector<unsigned long> t(2);
240  os << "Uvec\n";
241  t = DoubConv::dto2longs(defaultMean);
242  os << defaultMean << " " << t[0] << " " << t[1] << "\n";
243  t = DoubConv::dto2longs(defaultStdDev);
244  os << defaultStdDev << " " << t[0] << " " << t[1] << "\n";
245  if ( set ) {
246  t = DoubConv::dto2longs(nextGauss);
247  os << "nextGauss " << nextGauss << " " << t[0] << " " << t[1] << "\n";
248  } else {
249  os << "no_cached_nextGauss \n";
250  }
251  os.precision(prec);
252  return os;
253 } // put
254 
255 std::istream & RandGauss::get ( std::istream & is ) {
256  std::string inName;
257  is >> inName;
258  if (inName != name()) {
259  is.clear(std::ios::badbit | is.rdstate());
260  std::cerr << "Mismatch when expecting to read state of a "
261  << name() << " distribution\n"
262  << "Name found was " << inName
263  << "\nistream is left in the badbit state\n";
264  return is;
265  }
266  std::string c1;
267  std::string c2;
268  if (possibleKeywordInput(is, "Uvec", c1)) {
269  std::vector<unsigned long> t(2);
270  is >> defaultMean >> t[0] >> t[1]; defaultMean = DoubConv::longs2double(t);
271  is >> defaultStdDev>>t[0]>>t[1]; defaultStdDev = DoubConv::longs2double(t);
272  std::string ng;
273  is >> ng;
274  set = false;
275  if (ng == "nextGauss") {
276  is >> nextGauss >> t[0] >> t[1]; nextGauss = DoubConv::longs2double(t);
277  set = true;
278  }
279  return is;
280  }
281  // is >> c1 encompassed by possibleKeywordInput
282  is >> defaultMean >> c2 >> defaultStdDev;
283  if ( (!is) || (c1 != "Mean:") || (c2 != "Sigma:") ) {
284  std::cerr << "i/o problem while expecting to read state of a "
285  << name() << " distribution\n"
286  << "default mean and/or sigma could not be read\n";
287  return is;
288  }
289  is >> c1 >> c2 >> nextGauss;
290  if ( (!is) || (c1 != "RANDGAUSS") ) {
291  is.clear(std::ios::badbit | is.rdstate());
292  std::cerr << "Failure when reading caching state of RandGauss\n";
293  return is;
294  }
295  if (c2 == "CACHED_GAUSSIAN:") {
296  set = true;
297  } else if (c2 == "NO_CACHED_GAUSSIAN:") {
298  set = false;
299  } else {
300  is.clear(std::ios::badbit | is.rdstate());
301  std::cerr << "Unexpected caching state keyword of RandGauss:" << c2
302  << "\nistream is left in the badbit state\n";
303  }
304  return is;
305 } // get
306 
307  // Static save and restore to/from streams
308 
309 std::ostream & RandGauss::saveDistState ( std::ostream & os ) {
310  int prec = os.precision(20);
311  std::vector<unsigned long> t(2);
312  os << distributionName() << "\n";
313  os << "Uvec\n";
314  if ( getFlag() ) {
315  t = DoubConv::dto2longs(getVal());
316  os << "nextGauss_st " << getVal() << " " << t[0] << " " << t[1] << "\n";
317  } else {
318  os << "no_cached_nextGauss_st \n";
319  }
320  os.precision(prec);
321  return os;
322 }
323 
324 std::istream & RandGauss::restoreDistState ( std::istream & is ) {
325  std::string inName;
326  is >> inName;
327  if (inName != distributionName()) {
328  is.clear(std::ios::badbit | is.rdstate());
329  std::cerr << "Mismatch when expecting to read static state of a "
330  << distributionName() << " distribution\n"
331  << "Name found was " << inName
332  << "\nistream is left in the badbit state\n";
333  return is;
334  }
335  std::string c1;
336  std::string c2;
337  if (possibleKeywordInput(is, "Uvec", c1)) {
338  std::vector<unsigned long> t(2);
339  std::string ng;
340  is >> ng;
341  setFlag (false);
342  if (ng == "nextGauss_st") {
343  is >> nextGauss_st >> t[0] >> t[1];
344  nextGauss_st = DoubConv::longs2double(t);
345  setFlag (true);
346  }
347  return is;
348  }
349  // is >> c1 encompassed by possibleKeywordInput
350  is >> c2 >> nextGauss_st;
351  if ( (!is) || (c1 != "RANDGAUSS") ) {
352  is.clear(std::ios::badbit | is.rdstate());
353  std::cerr << "Failure when reading caching state of static RandGauss\n";
354  return is;
355  }
356  if (c2 == "CACHED_GAUSSIAN:") {
357  setFlag(true);
358  } else if (c2 == "NO_CACHED_GAUSSIAN:") {
359  setFlag(false);
360  } else {
361  is.clear(std::ios::badbit | is.rdstate());
362  std::cerr << "Unexpected caching state keyword of static RandGauss:" << c2
363  << "\nistream is left in the badbit state\n";
364  }
365  return is;
366 }
367 
368 std::ostream & RandGauss::saveFullState ( std::ostream & os ) {
369  HepRandom::saveFullState(os);
370  saveDistState(os);
371  return os;
372 }
373 
374 std::istream & RandGauss::restoreFullState ( std::istream & is ) {
375  HepRandom::restoreFullState(is);
376  restoreDistState(is);
377  return is;
378 }
379 
380 } // namespace CLHEP
381 
ThreeVector shoot(const G4int Ap, const G4int Af)
static c2_factory< G4double > c2
static const G4double fac
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:51
static const G4double c1