00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035 #include "CLHEP/Random/RandGauss.h"
00036 #include "CLHEP/Random/DoubConv.h"
00037 #include <string.h>
00038 #include <cmath>
00039
00040 namespace CLHEP {
00041
00042 std::string RandGauss::name() const {return "RandGauss";}
00043 HepRandomEngine & RandGauss::engine() {return *localEngine;}
00044
00045
00046 bool RandGauss::set_st = false;
00047 double RandGauss::nextGauss_st = 0.0;
00048
00049 RandGauss::~RandGauss() {
00050 }
00051
00052 double RandGauss::operator()() {
00053 return fire( defaultMean, defaultStdDev );
00054 }
00055
00056 double RandGauss::operator()( double mean, double stdDev ) {
00057 return fire( mean, stdDev );
00058 }
00059
00060 double RandGauss::shoot()
00061 {
00062
00063
00064
00065 if ( getFlag() ) {
00066 setFlag(false);
00067 double x = getVal();
00068 return x;
00069
00070 }
00071
00072 double r;
00073 double v1,v2,fac,val;
00074 HepRandomEngine* anEngine = HepRandom::getTheEngine();
00075
00076 do {
00077 v1 = 2.0 * anEngine->flat() - 1.0;
00078 v2 = 2.0 * anEngine->flat() - 1.0;
00079 r = v1*v1 + v2*v2;
00080 } while ( r > 1.0 );
00081
00082 fac = std::sqrt(-2.0*std::log(r)/r);
00083 val = v1*fac;
00084 setVal(val);
00085 setFlag(true);
00086 return v2*fac;
00087 }
00088
00089 void RandGauss::shootArray( const int size, double* vect,
00090 double mean, double stdDev )
00091 {
00092 for( double* v = vect; v != vect + size; ++v )
00093 *v = shoot(mean,stdDev);
00094 }
00095
00096 double RandGauss::shoot( HepRandomEngine* anEngine )
00097 {
00098
00099
00100
00101 if ( getFlag() ) {
00102 setFlag(false);
00103 return getVal();
00104 }
00105
00106 double r;
00107 double v1,v2,fac,val;
00108
00109 do {
00110 v1 = 2.0 * anEngine->flat() - 1.0;
00111 v2 = 2.0 * anEngine->flat() - 1.0;
00112 r = v1*v1 + v2*v2;
00113 } while ( r > 1.0 );
00114
00115 fac = std::sqrt( -2.0*std::log(r)/r);
00116 val = v1*fac;
00117 setVal(val);
00118 setFlag(true);
00119 return v2*fac;
00120 }
00121
00122 void RandGauss::shootArray( HepRandomEngine* anEngine,
00123 const int size, double* vect,
00124 double mean, double stdDev )
00125 {
00126 for( double* v = vect; v != vect + size; ++v )
00127 *v = shoot(anEngine,mean,stdDev);
00128 }
00129
00130 double RandGauss::normal()
00131 {
00132
00133
00134
00135 if ( set ) {
00136 set = false;
00137 return nextGauss;
00138 }
00139
00140 double r;
00141 double v1,v2,fac,val;
00142
00143 do {
00144 v1 = 2.0 * localEngine->flat() - 1.0;
00145 v2 = 2.0 * localEngine->flat() - 1.0;
00146 r = v1*v1 + v2*v2;
00147 } while ( r > 1.0 );
00148
00149 fac = std::sqrt(-2.0*std::log(r)/r);
00150 val = v1*fac;
00151 nextGauss = val;
00152 set = true;
00153 return v2*fac;
00154 }
00155
00156 void RandGauss::fireArray( const int size, double* vect)
00157 {
00158 for( double* v = vect; v != vect + size; ++v )
00159 *v = fire( defaultMean, defaultStdDev );
00160 }
00161
00162 void RandGauss::fireArray( const int size, double* vect,
00163 double mean, double stdDev )
00164 {
00165 for( double* v = vect; v != vect + size; ++v )
00166 *v = fire( mean, stdDev );
00167 }
00168
00169 void RandGauss::saveEngineStatus ( const char filename[] ) {
00170
00171
00172 getTheEngine()->saveStatus( filename );
00173
00174
00175
00176 std::ofstream outfile ( filename, std::ios::app );
00177
00178 if ( getFlag() ) {
00179 std::vector<unsigned long> t(2);
00180 t = DoubConv::dto2longs(getVal());
00181 outfile << "RANDGAUSS CACHED_GAUSSIAN: Uvec "
00182 << getVal() << " " << t[0] << " " << t[1] << "\n";
00183 } else {
00184 outfile << "RANDGAUSS NO_CACHED_GAUSSIAN: 0 \n" ;
00185 }
00186
00187 }
00188
00189 void RandGauss::restoreEngineStatus( const char filename[] ) {
00190
00191
00192 getTheEngine()->restoreStatus( filename );
00193
00194
00195
00196 std::ifstream infile ( filename, std::ios::in );
00197 if (!infile) return;
00198
00199 char inputword[] = "NO_KEYWORD ";
00200 while (true) {
00201 infile.width(13);
00202 infile >> inputword;
00203 if (strcmp(inputword,"RANDGAUSS")==0) break;
00204 if (infile.eof()) break;
00205
00206
00207
00208 }
00209
00210
00211
00212 if (strcmp(inputword,"RANDGAUSS")==0) {
00213 char setword[40];
00214 infile.width(39);
00215 infile >> setword;
00216 if (strcmp(setword,"CACHED_GAUSSIAN:") ==0) {
00217 if (possibleKeywordInput(infile, "Uvec", nextGauss_st)) {
00218 std::vector<unsigned long> t(2);
00219 infile >> nextGauss_st >> t[0] >> t[1];
00220 nextGauss_st = DoubConv::longs2double(t);
00221 }
00222
00223 setFlag(true);
00224 } else {
00225 setFlag(false);
00226 infile >> nextGauss_st;
00227 }
00228 } else {
00229 setFlag(false);
00230 }
00231
00232 }
00233
00234
00235
00236 std::ostream & RandGauss::put ( std::ostream & os ) const {
00237 os << name() << "\n";
00238 int prec = os.precision(20);
00239 std::vector<unsigned long> t(2);
00240 os << "Uvec\n";
00241 t = DoubConv::dto2longs(defaultMean);
00242 os << defaultMean << " " << t[0] << " " << t[1] << "\n";
00243 t = DoubConv::dto2longs(defaultStdDev);
00244 os << defaultStdDev << " " << t[0] << " " << t[1] << "\n";
00245 if ( set ) {
00246 t = DoubConv::dto2longs(nextGauss);
00247 os << "nextGauss " << nextGauss << " " << t[0] << " " << t[1] << "\n";
00248 } else {
00249 os << "no_cached_nextGauss \n";
00250 }
00251 os.precision(prec);
00252 return os;
00253 }
00254
00255 std::istream & RandGauss::get ( std::istream & is ) {
00256 std::string inName;
00257 is >> inName;
00258 if (inName != name()) {
00259 is.clear(std::ios::badbit | is.rdstate());
00260 std::cerr << "Mismatch when expecting to read state of a "
00261 << name() << " distribution\n"
00262 << "Name found was " << inName
00263 << "\nistream is left in the badbit state\n";
00264 return is;
00265 }
00266 std::string c1;
00267 std::string c2;
00268 if (possibleKeywordInput(is, "Uvec", c1)) {
00269 std::vector<unsigned long> t(2);
00270 is >> defaultMean >> t[0] >> t[1]; defaultMean = DoubConv::longs2double(t);
00271 is >> defaultStdDev>>t[0]>>t[1]; defaultStdDev = DoubConv::longs2double(t);
00272 std::string ng;
00273 is >> ng;
00274 set = false;
00275 if (ng == "nextGauss") {
00276 is >> nextGauss >> t[0] >> t[1]; nextGauss = DoubConv::longs2double(t);
00277 set = true;
00278 }
00279 return is;
00280 }
00281
00282 is >> defaultMean >> c2 >> defaultStdDev;
00283 if ( (!is) || (c1 != "Mean:") || (c2 != "Sigma:") ) {
00284 std::cerr << "i/o problem while expecting to read state of a "
00285 << name() << " distribution\n"
00286 << "default mean and/or sigma could not be read\n";
00287 return is;
00288 }
00289 is >> c1 >> c2 >> nextGauss;
00290 if ( (!is) || (c1 != "RANDGAUSS") ) {
00291 is.clear(std::ios::badbit | is.rdstate());
00292 std::cerr << "Failure when reading caching state of RandGauss\n";
00293 return is;
00294 }
00295 if (c2 == "CACHED_GAUSSIAN:") {
00296 set = true;
00297 } else if (c2 == "NO_CACHED_GAUSSIAN:") {
00298 set = false;
00299 } else {
00300 is.clear(std::ios::badbit | is.rdstate());
00301 std::cerr << "Unexpected caching state keyword of RandGauss:" << c2
00302 << "\nistream is left in the badbit state\n";
00303 }
00304 return is;
00305 }
00306
00307
00308
00309 std::ostream & RandGauss::saveDistState ( std::ostream & os ) {
00310 int prec = os.precision(20);
00311 std::vector<unsigned long> t(2);
00312 os << distributionName() << "\n";
00313 os << "Uvec\n";
00314 if ( getFlag() ) {
00315 t = DoubConv::dto2longs(getVal());
00316 os << "nextGauss_st " << getVal() << " " << t[0] << " " << t[1] << "\n";
00317 } else {
00318 os << "no_cached_nextGauss_st \n";
00319 }
00320 os.precision(prec);
00321 return os;
00322 }
00323
00324 std::istream & RandGauss::restoreDistState ( std::istream & is ) {
00325 std::string inName;
00326 is >> inName;
00327 if (inName != distributionName()) {
00328 is.clear(std::ios::badbit | is.rdstate());
00329 std::cerr << "Mismatch when expecting to read static state of a "
00330 << distributionName() << " distribution\n"
00331 << "Name found was " << inName
00332 << "\nistream is left in the badbit state\n";
00333 return is;
00334 }
00335 std::string c1;
00336 std::string c2;
00337 if (possibleKeywordInput(is, "Uvec", c1)) {
00338 std::vector<unsigned long> t(2);
00339 std::string ng;
00340 is >> ng;
00341 setFlag (false);
00342 if (ng == "nextGauss_st") {
00343 is >> nextGauss_st >> t[0] >> t[1];
00344 nextGauss_st = DoubConv::longs2double(t);
00345 setFlag (true);
00346 }
00347 return is;
00348 }
00349
00350 is >> c2 >> nextGauss_st;
00351 if ( (!is) || (c1 != "RANDGAUSS") ) {
00352 is.clear(std::ios::badbit | is.rdstate());
00353 std::cerr << "Failure when reading caching state of static RandGauss\n";
00354 return is;
00355 }
00356 if (c2 == "CACHED_GAUSSIAN:") {
00357 setFlag(true);
00358 } else if (c2 == "NO_CACHED_GAUSSIAN:") {
00359 setFlag(false);
00360 } else {
00361 is.clear(std::ios::badbit | is.rdstate());
00362 std::cerr << "Unexpected caching state keyword of static RandGauss:" << c2
00363 << "\nistream is left in the badbit state\n";
00364 }
00365 return is;
00366 }
00367
00368 std::ostream & RandGauss::saveFullState ( std::ostream & os ) {
00369 HepRandom::saveFullState(os);
00370 saveDistState(os);
00371 return os;
00372 }
00373
00374 std::istream & RandGauss::restoreFullState ( std::istream & is ) {
00375 HepRandom::restoreFullState(is);
00376 restoreDistState(is);
00377 return is;
00378 }
00379
00380 }
00381