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