CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: RandGauss.cc,v 1.6 2010/06/16 17:24: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> // for strcmp 00039 #include <cmath> // for std::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 } 00052 00053 double RandGauss::operator()() { 00054 return fire( defaultMean, defaultStdDev ); 00055 } 00056 00057 double RandGauss::operator()( double mean, double stdDev ) { 00058 return fire( mean, stdDev ); 00059 } 00060 00061 double RandGauss::shoot() 00062 { 00063 // Gaussian random numbers are generated two at the time, so every other 00064 // time this is called we just return a number generated the time before. 00065 00066 if ( getFlag() ) { 00067 setFlag(false); 00068 double x = getVal(); 00069 return x; 00070 // return getVal(); 00071 } 00072 00073 double r; 00074 double v1,v2,fac,val; 00075 HepRandomEngine* anEngine = HepRandom::getTheEngine(); 00076 00077 do { 00078 v1 = 2.0 * anEngine->flat() - 1.0; 00079 v2 = 2.0 * anEngine->flat() - 1.0; 00080 r = v1*v1 + v2*v2; 00081 } while ( r > 1.0 ); 00082 00083 fac = std::sqrt(-2.0*std::log(r)/r); 00084 val = v1*fac; 00085 setVal(val); 00086 setFlag(true); 00087 return v2*fac; 00088 } 00089 00090 void RandGauss::shootArray( const int size, double* vect, 00091 double mean, double stdDev ) 00092 { 00093 for( double* v = vect; v != vect + size; ++v ) 00094 *v = shoot(mean,stdDev); 00095 } 00096 00097 double RandGauss::shoot( HepRandomEngine* anEngine ) 00098 { 00099 // Gaussian random numbers are generated two at the time, so every other 00100 // time this is called we just return a number generated the time before. 00101 00102 if ( getFlag() ) { 00103 setFlag(false); 00104 return getVal(); 00105 } 00106 00107 double r; 00108 double v1,v2,fac,val; 00109 00110 do { 00111 v1 = 2.0 * anEngine->flat() - 1.0; 00112 v2 = 2.0 * anEngine->flat() - 1.0; 00113 r = v1*v1 + v2*v2; 00114 } while ( r > 1.0 ); 00115 00116 fac = std::sqrt( -2.0*std::log(r)/r); 00117 val = v1*fac; 00118 setVal(val); 00119 setFlag(true); 00120 return v2*fac; 00121 } 00122 00123 void RandGauss::shootArray( HepRandomEngine* anEngine, 00124 const int size, double* vect, 00125 double mean, double stdDev ) 00126 { 00127 for( double* v = vect; v != vect + size; ++v ) 00128 *v = shoot(anEngine,mean,stdDev); 00129 } 00130 00131 double RandGauss::normal() 00132 { 00133 // Gaussian random numbers are generated two at the time, so every other 00134 // time this is called we just return a number generated the time before. 00135 00136 if ( set ) { 00137 set = false; 00138 return nextGauss; 00139 } 00140 00141 double r; 00142 double v1,v2,fac,val; 00143 00144 do { 00145 v1 = 2.0 * localEngine->flat() - 1.0; 00146 v2 = 2.0 * localEngine->flat() - 1.0; 00147 r = v1*v1 + v2*v2; 00148 } while ( r > 1.0 ); 00149 00150 fac = std::sqrt(-2.0*std::log(r)/r); 00151 val = v1*fac; 00152 nextGauss = val; 00153 set = true; 00154 return v2*fac; 00155 } 00156 00157 void RandGauss::fireArray( const int size, double* vect) 00158 { 00159 for( double* v = vect; v != vect + size; ++v ) 00160 *v = fire( defaultMean, defaultStdDev ); 00161 } 00162 00163 void RandGauss::fireArray( const int size, double* vect, 00164 double mean, double stdDev ) 00165 { 00166 for( double* v = vect; v != vect + size; ++v ) 00167 *v = fire( mean, stdDev ); 00168 } 00169 00170 void RandGauss::saveEngineStatus ( const char filename[] ) { 00171 00172 // First save the engine status just like the base class would do: 00173 getTheEngine()->saveStatus( filename ); 00174 00175 // Now append the cached variate, if any: 00176 00177 std::ofstream outfile ( filename, std::ios::app ); 00178 00179 if ( getFlag() ) { 00180 std::vector<unsigned long> t(2); 00181 t = DoubConv::dto2longs(getVal()); 00182 outfile << "RANDGAUSS CACHED_GAUSSIAN: Uvec " 00183 << getVal() << " " << t[0] << " " << t[1] << "\n"; 00184 } else { 00185 outfile << "RANDGAUSS NO_CACHED_GAUSSIAN: 0 \n" ; 00186 } 00187 00188 } // saveEngineStatus 00189 00190 void RandGauss::restoreEngineStatus( const char filename[] ) { 00191 00192 // First restore the engine status just like the base class would do: 00193 getTheEngine()->restoreStatus( filename ); 00194 00195 // Now find the line describing the cached variate: 00196 00197 std::ifstream infile ( filename, std::ios::in ); 00198 if (!infile) return; 00199 00200 char inputword[] = "NO_KEYWORD "; // leaves room for 14 characters plus \0 00201 while (true) { 00202 infile.width(13); 00203 infile >> inputword; 00204 if (strcmp(inputword,"RANDGAUSS")==0) break; 00205 if (infile.eof()) break; 00206 // If the file ends without the RANDGAUSS line, that means this 00207 // was a file produced by an earlier version of RandGauss. We will 00208 // replicated the old behavior in that case: set_st is cleared. 00209 } 00210 00211 // Then read and use the caching info: 00212 00213 if (strcmp(inputword,"RANDGAUSS")==0) { 00214 char setword[40]; // the longest, staticFirstUnusedBit: has length 21 00215 infile.width(39); 00216 infile >> setword; // setword should be CACHED_GAUSSIAN: 00217 if (strcmp(setword,"CACHED_GAUSSIAN:") ==0) { 00218 if (possibleKeywordInput(infile, "Uvec", nextGauss_st)) { 00219 std::vector<unsigned long> t(2); 00220 infile >> nextGauss_st >> t[0] >> t[1]; 00221 nextGauss_st = DoubConv::longs2double(t); 00222 } 00223 // is >> nextGauss_st encompassed by possibleKeywordInput 00224 setFlag(true); 00225 } else { 00226 setFlag(false); 00227 infile >> nextGauss_st; // because a 0 will have been output 00228 } 00229 } else { 00230 setFlag(false); 00231 } 00232 00233 } // restoreEngineStatus 00234 00235 // Save and restore to/from streams 00236 00237 std::ostream & RandGauss::put ( std::ostream & os ) const { 00238 os << name() << "\n"; 00239 int prec = os.precision(20); 00240 std::vector<unsigned long> t(2); 00241 os << "Uvec\n"; 00242 t = DoubConv::dto2longs(defaultMean); 00243 os << defaultMean << " " << t[0] << " " << t[1] << "\n"; 00244 t = DoubConv::dto2longs(defaultStdDev); 00245 os << defaultStdDev << " " << t[0] << " " << t[1] << "\n"; 00246 if ( set ) { 00247 t = DoubConv::dto2longs(nextGauss); 00248 os << "nextGauss " << nextGauss << " " << t[0] << " " << t[1] << "\n"; 00249 #ifdef TRACE_IO 00250 std::cout << "put(): nextGauss = " << nextGauss << "\n"; 00251 #endif 00252 } else { 00253 #ifdef TRACE_IO 00254 std::cout << "put(): No nextGauss \n"; 00255 #endif 00256 os << "no_cached_nextGauss \n"; 00257 } 00258 os.precision(prec); 00259 return os; 00260 } // put 00261 00262 std::istream & RandGauss::get ( std::istream & is ) { 00263 std::string inName; 00264 is >> inName; 00265 if (inName != name()) { 00266 is.clear(std::ios::badbit | is.rdstate()); 00267 std::cerr << "Mismatch when expecting to read state of a " 00268 << name() << " distribution\n" 00269 << "Name found was " << inName 00270 << "\nistream is left in the badbit state\n"; 00271 return is; 00272 } 00273 std::string c1; 00274 std::string c2; 00275 if (possibleKeywordInput(is, "Uvec", c1)) { 00276 std::vector<unsigned long> t(2); 00277 is >> defaultMean >> t[0] >> t[1]; defaultMean = DoubConv::longs2double(t); 00278 is >> defaultStdDev>>t[0]>>t[1]; defaultStdDev = DoubConv::longs2double(t); 00279 std::string ng; 00280 is >> ng; 00281 set = false; 00282 #ifdef TRACE_IO 00283 if (ng != "nextGauss") 00284 std::cout << "get(): ng = " << ng << "\n"; 00285 #endif 00286 if (ng == "nextGauss") { 00287 is >> nextGauss >> t[0] >> t[1]; nextGauss = DoubConv::longs2double(t); 00288 #ifdef TRACE_IO 00289 std::cout << "get(): nextGauss read back as " << nextGauss << "\n"; 00290 #endif 00291 set = true; 00292 } 00293 return is; 00294 } 00295 // is >> c1 encompassed by possibleKeywordInput 00296 is >> defaultMean >> c2 >> defaultStdDev; 00297 if ( (!is) || (c1 != "Mean:") || (c2 != "Sigma:") ) { 00298 std::cerr << "i/o problem while expecting to read state of a " 00299 << name() << " distribution\n" 00300 << "default mean and/or sigma could not be read\n"; 00301 return is; 00302 } 00303 is >> c1 >> c2 >> nextGauss; 00304 if ( (!is) || (c1 != "RANDGAUSS") ) { 00305 is.clear(std::ios::badbit | is.rdstate()); 00306 std::cerr << "Failure when reading caching state of RandGauss\n"; 00307 return is; 00308 } 00309 if (c2 == "CACHED_GAUSSIAN:") { 00310 set = true; 00311 } else if (c2 == "NO_CACHED_GAUSSIAN:") { 00312 set = false; 00313 } else { 00314 is.clear(std::ios::badbit | is.rdstate()); 00315 std::cerr << "Unexpected caching state keyword of RandGauss:" << c2 00316 << "\nistream is left in the badbit state\n"; 00317 } 00318 return is; 00319 } // get 00320 00321 // Static save and restore to/from streams 00322 00323 std::ostream & RandGauss::saveDistState ( std::ostream & os ) { 00324 int prec = os.precision(20); 00325 std::vector<unsigned long> t(2); 00326 os << distributionName() << "\n"; 00327 os << "Uvec\n"; 00328 if ( getFlag() ) { 00329 t = DoubConv::dto2longs(getVal()); 00330 os << "nextGauss_st " << getVal() << " " << t[0] << " " << t[1] << "\n"; 00331 } else { 00332 os << "no_cached_nextGauss_st \n"; 00333 } 00334 os.precision(prec); 00335 return os; 00336 } 00337 00338 std::istream & RandGauss::restoreDistState ( std::istream & is ) { 00339 std::string inName; 00340 is >> inName; 00341 if (inName != distributionName()) { 00342 is.clear(std::ios::badbit | is.rdstate()); 00343 std::cerr << "Mismatch when expecting to read static state of a " 00344 << distributionName() << " distribution\n" 00345 << "Name found was " << inName 00346 << "\nistream is left in the badbit state\n"; 00347 return is; 00348 } 00349 std::string c1; 00350 std::string c2; 00351 if (possibleKeywordInput(is, "Uvec", c1)) { 00352 std::vector<unsigned long> t(2); 00353 std::string ng; 00354 is >> ng; 00355 setFlag (false); 00356 if (ng == "nextGauss_st") { 00357 is >> nextGauss_st >> t[0] >> t[1]; 00358 nextGauss_st = DoubConv::longs2double(t); 00359 setFlag (true); 00360 } 00361 return is; 00362 } 00363 // is >> c1 encompassed by possibleKeywordInput 00364 is >> c2 >> nextGauss_st; 00365 if ( (!is) || (c1 != "RANDGAUSS") ) { 00366 is.clear(std::ios::badbit | is.rdstate()); 00367 std::cerr << "Failure when reading caching state of static RandGauss\n"; 00368 return is; 00369 } 00370 if (c2 == "CACHED_GAUSSIAN:") { 00371 setFlag(true); 00372 } else if (c2 == "NO_CACHED_GAUSSIAN:") { 00373 setFlag(false); 00374 } else { 00375 is.clear(std::ios::badbit | is.rdstate()); 00376 std::cerr << "Unexpected caching state keyword of static RandGauss:" << c2 00377 << "\nistream is left in the badbit state\n"; 00378 } 00379 return is; 00380 } 00381 00382 std::ostream & RandGauss::saveFullState ( std::ostream & os ) { 00383 HepRandom::saveFullState(os); 00384 saveDistState(os); 00385 return os; 00386 } 00387 00388 std::istream & RandGauss::restoreFullState ( std::istream & is ) { 00389 HepRandom::restoreFullState(is); 00390 restoreDistState(is); 00391 return is; 00392 } 00393 00394 } // namespace CLHEP 00395