CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: RandGauss.cc,v 1.4.4.2 2005/04/15 16:32:53 garren Exp $ 00002 // -*- C++ -*- 00003 // 00004 // ----------------------------------------------------------------------- 00005 // HEP Random 00006 // --- RandGauss --- 00007 // class implementation file 00008 // ----------------------------------------------------------------------- 00009 // This file is part of Geant4 (simulation toolkit for HEP). 00010 00011 // ======================================================================= 00012 // Gabriele Cosmo - Created: 5th September 1995 00013 // - Added methods to shoot arrays: 28th July 1997 00014 // J.Marraffino - Added default arguments as attributes and 00015 // operator() with arguments. Introduced method normal() 00016 // for computation in fire(): 16th Feb 1998 00017 // Gabriele Cosmo - Relocated static data from HepRandom: 5th Jan 1999 00018 // M Fischler - Copy constructor should supply right engine to HepRandom: 00019 // 1/26/00. 00020 // M Fischler - Workaround for problem of non-reproducing saveEngineStatus 00021 // by saving cached gaussian. March 2000. 00022 // M Fischler - Avoiding hang when file not found in restoreEngineStatus 00023 // 12/3/04 00024 // M Fischler - put and get to/from streams 12/8/04 00025 // M Fischler - save and restore dist to streams 12/20/04 00026 // M Fischler - put/get to/from streams uses pairs of ulongs when 00027 // storing doubles avoid problems with precision. 00028 // Similarly for saveEngineStatus and RestoreEngineStatus 00029 // and for save/restore distState 00030 // Care was taken that old-form output can still be read back. 00031 // 4/14/05 00032 // 00033 // ======================================================================= 00034 00035 #include "CLHEP/Random/defs.h" 00036 #include "CLHEP/Random/RandGauss.h" 00037 #include "CLHEP/Random/DoubConv.hh" 00038 #include <string.h> 00039 #include <cmath> // for log() 00040 00041 namespace CLHEP { 00042 00043 std::string RandGauss::name() const {return "RandGauss";} 00044 HepRandomEngine & RandGauss::engine() {return *localEngine;} 00045 00046 // Initialisation of static data 00047 bool RandGauss::set_st = false; 00048 double RandGauss::nextGauss_st = 0.0; 00049 00050 RandGauss::~RandGauss() { 00051 if ( deleteEngine ) delete localEngine; 00052 } 00053 00054 RandGauss::RandGauss(const RandGauss& right) 00055 : HepRandom(right.getTheEngine()), 00056 defaultMean(right.defaultMean), defaultStdDev(right.defaultStdDev) 00057 {;} 00058 00059 double RandGauss::operator()() { 00060 return fire( defaultMean, defaultStdDev ); 00061 } 00062 00063 double RandGauss::operator()( double mean, double stdDev ) { 00064 return fire( mean, stdDev ); 00065 } 00066 00067 double RandGauss::shoot() 00068 { 00069 // Gaussian random numbers are generated two at the time, so every other 00070 // time this is called we just return a number generated the time before. 00071 00072 if ( getFlag() ) { 00073 setFlag(false); 00074 double x = getVal(); 00075 return x; 00076 // return getVal(); 00077 } 00078 00079 double r; 00080 double v1,v2,fac,val; 00081 HepRandomEngine* anEngine = HepRandom::getTheEngine(); 00082 00083 do { 00084 v1 = 2.0 * anEngine->flat() - 1.0; 00085 v2 = 2.0 * anEngine->flat() - 1.0; 00086 r = v1*v1 + v2*v2; 00087 } while ( r > 1.0 ); 00088 00089 fac = sqrt(-2.0*log(r)/r); 00090 val = v1*fac; 00091 setVal(val); 00092 setFlag(true); 00093 return v2*fac; 00094 } 00095 00096 void RandGauss::shootArray( const int size, double* vect, 00097 double mean, double stdDev ) 00098 { 00099 int i; 00100 00101 for (i=0; i<size; ++i) 00102 vect[i] = shoot(mean,stdDev); 00103 } 00104 00105 double RandGauss::shoot( HepRandomEngine* anEngine ) 00106 { 00107 // Gaussian random numbers are generated two at the time, so every other 00108 // time this is called we just return a number generated the time before. 00109 00110 if ( getFlag() ) { 00111 setFlag(false); 00112 return getVal(); 00113 } 00114 00115 double r; 00116 double v1,v2,fac,val; 00117 00118 do { 00119 v1 = 2.0 * anEngine->flat() - 1.0; 00120 v2 = 2.0 * anEngine->flat() - 1.0; 00121 r = v1*v1 + v2*v2; 00122 } while ( r > 1.0 ); 00123 00124 fac = sqrt( -2.0*log(r)/r); 00125 val = v1*fac; 00126 setVal(val); 00127 setFlag(true); 00128 return v2*fac; 00129 } 00130 00131 void RandGauss::shootArray( HepRandomEngine* anEngine, 00132 const int size, double* vect, 00133 double mean, double stdDev ) 00134 { 00135 int i; 00136 00137 for (i=0; i<size; ++i) 00138 vect[i] = shoot(anEngine,mean,stdDev); 00139 } 00140 00141 double RandGauss::normal() 00142 { 00143 // Gaussian random numbers are generated two at the time, so every other 00144 // time this is called we just return a number generated the time before. 00145 00146 if ( set ) { 00147 set = false; 00148 return nextGauss; 00149 } 00150 00151 double r; 00152 double v1,v2,fac,val; 00153 00154 do { 00155 v1 = 2.0 * localEngine->flat() - 1.0; 00156 v2 = 2.0 * localEngine->flat() - 1.0; 00157 r = v1*v1 + v2*v2; 00158 } while ( r > 1.0 ); 00159 00160 fac = sqrt(-2.0*log(r)/r); 00161 val = v1*fac; 00162 nextGauss = val; 00163 set = true; 00164 return v2*fac; 00165 } 00166 00167 void RandGauss::fireArray( const int size, double* vect) 00168 { 00169 int i; 00170 00171 for (i=0; i<size; ++i) 00172 vect[i] = fire( defaultMean, defaultStdDev ); 00173 } 00174 00175 void RandGauss::fireArray( const int size, double* vect, 00176 double mean, double stdDev ) 00177 { 00178 int i; 00179 00180 for (i=0; i<size; ++i) 00181 vect[i] = fire( mean, stdDev ); 00182 } 00183 00184 void RandGauss::saveEngineStatus ( const char filename[] ) { 00185 00186 // First save the engine status just like the base class would do: 00187 getTheEngine()->saveStatus( filename ); 00188 00189 // Now append the cached variate, if any: 00190 00191 std::ofstream outfile ( filename, std::ios::app ); 00192 00193 if ( getFlag() ) { 00194 std::vector<unsigned long> t(2); 00195 t = DoubConv::dto2longs(getVal()); 00196 outfile << "RANDGAUSS CACHED_GAUSSIAN: Uvec " 00197 << getVal() << " " << t[0] << " " << t[1] << "\n"; 00198 } else { 00199 outfile << "RANDGAUSS NO_CACHED_GAUSSIAN: 0 \n" ; 00200 } 00201 00202 } // saveEngineStatus 00203 00204 void RandGauss::restoreEngineStatus( const char filename[] ) { 00205 00206 // First restore the engine status just like the base class would do: 00207 getTheEngine()->restoreStatus( filename ); 00208 00209 // Now find the line describing the cached variate: 00210 00211 std::ifstream infile ( filename, std::ios::in ); 00212 if (!infile) return; 00213 00214 char inputword[] = "NO_KEYWORD "; // leaves room for 14 characters plus \0 00215 while (true) { 00216 infile.width(13); 00217 infile >> inputword; 00218 if (strcmp(inputword,"RANDGAUSS")==0) break; 00219 if (infile.eof()) break; 00220 // If the file ends without the RANDGAUSS line, that means this 00221 // was a file produced by an earlier version of RandGauss. We will 00222 // replicated the old behavior in that case: set_st is cleared. 00223 } 00224 00225 // Then read and use the caching info: 00226 00227 if (strcmp(inputword,"RANDGAUSS")==0) { 00228 char setword[40]; // the longest, staticFirstUnusedBit: has length 21 00229 infile.width(39); 00230 infile >> setword; // setword should be CACHED_GAUSSIAN: 00231 if (strcmp(setword,"CACHED_GAUSSIAN:") ==0) { 00232 if (possibleKeywordInput(infile, "Uvec", nextGauss_st)) { 00233 std::vector<unsigned long> t(2); 00234 infile >> nextGauss_st >> t[0] >> t[1]; 00235 nextGauss_st = DoubConv::longs2double(t); 00236 } 00237 // is >> nextGauss_st encompassed by possibleKeywordInput 00238 setFlag(true); 00239 } else { 00240 setFlag(false); 00241 infile >> nextGauss_st; // because a 0 will have been output 00242 } 00243 } else { 00244 setFlag(false); 00245 } 00246 00247 } // restoreEngineStatus 00248 00249 // Save and restore to/from streams 00250 00251 std::ostream & RandGauss::put ( std::ostream & os ) const { 00252 os << name() << "\n"; 00253 int prec = os.precision(20); 00254 std::vector<unsigned long> t(2); 00255 os << "Uvec\n"; 00256 t = DoubConv::dto2longs(defaultMean); 00257 os << defaultMean << " " << t[0] << " " << t[1] << "\n"; 00258 t = DoubConv::dto2longs(defaultStdDev); 00259 os << defaultStdDev << " " << t[0] << " " << t[1] << "\n"; 00260 if ( set ) { 00261 t = DoubConv::dto2longs(nextGauss); 00262 os << "nextGauss " << nextGauss << " " << t[0] << " " << t[1] << "\n"; 00263 #ifdef TRACE_IO 00264 std::cout << "put(): nextGauss = " << nextGauss << "\n"; 00265 #endif 00266 } else { 00267 #ifdef TRACE_IO 00268 std::cout << "put(): No nextGauss \n"; 00269 #endif 00270 os << "no_cached_nextGauss \n"; 00271 } 00272 os.precision(prec); 00273 return os; 00274 } // put 00275 00276 std::istream & RandGauss::get ( std::istream & is ) { 00277 std::string inName; 00278 is >> inName; 00279 if (inName != name()) { 00280 is.clear(std::ios::badbit | is.rdstate()); 00281 std::cerr << "Mismatch when expecting to read state of a " 00282 << name() << " distribution\n" 00283 << "Name found was " << inName 00284 << "\nistream is left in the badbit state\n"; 00285 return is; 00286 } 00287 std::string c1; 00288 std::string c2; 00289 if (possibleKeywordInput(is, "Uvec", c1)) { 00290 std::vector<unsigned long> t(2); 00291 is >> defaultMean >> t[0] >> t[1]; defaultMean = DoubConv::longs2double(t); 00292 is >> defaultStdDev>>t[0]>>t[1]; defaultStdDev = DoubConv::longs2double(t); 00293 std::string ng; 00294 is >> ng; 00295 set = false; 00296 #ifdef TRACE_IO 00297 if (ng != "nextGauss") 00298 std::cout << "get(): ng = " << ng << "\n"; 00299 #endif 00300 if (ng == "nextGauss") { 00301 is >> nextGauss >> t[0] >> t[1]; nextGauss = DoubConv::longs2double(t); 00302 #ifdef TRACE_IO 00303 std::cout << "get(): nextGauss read back as " << nextGauss << "\n"; 00304 #endif 00305 set = true; 00306 } 00307 return is; 00308 } 00309 // is >> c1 encompassed by possibleKeywordInput 00310 is >> defaultMean >> c2 >> defaultStdDev; 00311 if ( (!is) || (c1 != "Mean:") || (c2 != "Sigma:") ) { 00312 std::cerr << "i/o problem while expecting to read state of a " 00313 << name() << " distribution\n" 00314 << "default mean and/or sigma could not be read\n"; 00315 return is; 00316 } 00317 is >> c1 >> c2 >> nextGauss; 00318 if ( (!is) || (c1 != "RANDGAUSS") ) { 00319 is.clear(std::ios::badbit | is.rdstate()); 00320 std::cerr << "Failure when reading caching state of RandGauss\n"; 00321 return is; 00322 } 00323 if (c2 == "CACHED_GAUSSIAN:") { 00324 set = true; 00325 } else if (c2 == "NO_CACHED_GAUSSIAN:") { 00326 set = false; 00327 } else { 00328 is.clear(std::ios::badbit | is.rdstate()); 00329 std::cerr << "Unexpected caching state keyword of RandGauss:" << c2 00330 << "\nistream is left in the badbit state\n"; 00331 } 00332 return is; 00333 } // get 00334 00335 // Static save and restore to/from streams 00336 00337 std::ostream & RandGauss::saveDistState ( std::ostream & os ) { 00338 int prec = os.precision(20); 00339 std::vector<unsigned long> t(2); 00340 os << distributionName() << "\n"; 00341 os << "Uvec\n"; 00342 if ( getFlag() ) { 00343 t = DoubConv::dto2longs(getVal()); 00344 os << "nextGauss_st " << getVal() << " " << t[0] << " " << t[1] << "\n"; 00345 } else { 00346 os << "no_cached_nextGauss_st \n"; 00347 } 00348 os.precision(prec); 00349 return os; 00350 } 00351 00352 std::istream & RandGauss::restoreDistState ( std::istream & is ) { 00353 std::string inName; 00354 is >> inName; 00355 if (inName != distributionName()) { 00356 is.clear(std::ios::badbit | is.rdstate()); 00357 std::cerr << "Mismatch when expecting to read static state of a " 00358 << distributionName() << " distribution\n" 00359 << "Name found was " << inName 00360 << "\nistream is left in the badbit state\n"; 00361 return is; 00362 } 00363 std::string c1; 00364 std::string c2; 00365 if (possibleKeywordInput(is, "Uvec", c1)) { 00366 std::vector<unsigned long> t(2); 00367 std::string ng; 00368 is >> ng; 00369 setFlag (false); 00370 if (ng == "nextGauss_st") { 00371 is >> nextGauss_st >> t[0] >> t[1]; 00372 nextGauss_st = DoubConv::longs2double(t); 00373 setFlag (true); 00374 } 00375 return is; 00376 } 00377 // is >> c1 encompassed by possibleKeywordInput 00378 is >> c2 >> nextGauss_st; 00379 if ( (!is) || (c1 != "RANDGAUSS") ) { 00380 is.clear(std::ios::badbit | is.rdstate()); 00381 std::cerr << "Failure when reading caching state of static RandGauss\n"; 00382 return is; 00383 } 00384 if (c2 == "CACHED_GAUSSIAN:") { 00385 setFlag(true); 00386 } else if (c2 == "NO_CACHED_GAUSSIAN:") { 00387 setFlag(false); 00388 } else { 00389 is.clear(std::ios::badbit | is.rdstate()); 00390 std::cerr << "Unexpected caching state keyword of static RandGauss:" << c2 00391 << "\nistream is left in the badbit state\n"; 00392 } 00393 return is; 00394 } 00395 00396 std::ostream & RandGauss::saveFullState ( std::ostream & os ) { 00397 HepRandom::saveFullState(os); 00398 saveDistState(os); 00399 return os; 00400 } 00401 00402 std::istream & RandGauss::restoreFullState ( std::istream & is ) { 00403 HepRandom::restoreFullState(is); 00404 restoreDistState(is); 00405 return is; 00406 } 00407 00408 } // namespace CLHEP 00409