CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: RandStudentT.cc,v 1.4.4.2 2005/04/15 16:32:53 garren Exp $ 00002 // -*- C++ -*- 00003 // 00004 // ----------------------------------------------------------------------- 00005 // HEP Random 00006 // --- RandStudentT --- 00007 // class implementation file 00008 // ----------------------------------------------------------------------- 00009 00010 // ======================================================================= 00011 // John Marraffino - Created: 12th May 1998 00012 // G.Cosmo - Fixed minor bug on inline definition for shoot() 00013 // methods : 20th Aug 1998 00014 // M Fischler - put and get to/from streams 12/13/04 00015 // M Fischler - put/get to/from streams uses pairs of ulongs when 00016 // + storing doubles avoid problems with precision 00017 // 4/14/05 00018 // ======================================================================= 00019 00020 #include <float.h> 00021 #include <cmath> // for log() exp() 00022 #include "CLHEP/Random/defs.h" 00023 #include "CLHEP/Random/RandStudentT.h" 00024 #include "CLHEP/Random/DoubConv.hh" 00025 00026 namespace CLHEP { 00027 00028 std::string RandStudentT::name() const {return "RandStudentT";} 00029 HepRandomEngine & RandStudentT::engine() {return *localEngine;} 00030 00031 RandStudentT::~RandStudentT() { 00032 if ( deleteEngine ) delete localEngine; 00033 } 00034 00035 RandStudentT::RandStudentT(const RandStudentT& right) 00036 : defaultA(right.defaultA) 00037 {;} 00038 00039 double RandStudentT::operator()() { 00040 return fire( defaultA ); 00041 } 00042 00043 double RandStudentT::operator()( double a ) { 00044 return fire( a ); 00045 } 00046 00047 double RandStudentT::shoot( double a ) { 00048 /****************************************************************** 00049 * * 00050 * Student-t Distribution - Polar Method * 00051 * * 00052 ****************************************************************** 00053 * * 00054 * The polar method of Box/Muller for generating Normal variates * 00055 * is adapted to the Student-t distribution. The two generated * 00056 * variates are not independent and the expected no. of uniforms * 00057 * per variate is 2.5464. * 00058 * * 00059 ****************************************************************** 00060 * * 00061 * FUNCTION : - tpol samples a random number from the * 00062 * Student-t distribution with parameter a > 0. * 00063 * REFERENCE : - R.W. Bailey (1994): Polar generation of random * 00064 * variates with the t-distribution, Mathematics * 00065 * of Computation 62, 779-781. * 00066 * SUBPROGRAM : - ... (0,1)-Uniform generator * 00067 * * 00068 * * 00069 * Implemented by F. Niederl, 1994 * 00070 ******************************************************************/ 00071 00072 double u,v,w; 00073 00074 // Check for valid input value 00075 00076 if( a < 0.0) return (DBL_MAX); 00077 00078 do 00079 { 00080 u = 2.0 * HepRandom::getTheEngine()->flat() - 1.0; 00081 v = 2.0 * HepRandom::getTheEngine()->flat() - 1.0; 00082 } 00083 while ((w = u * u + v * v) > 1.0); 00084 00085 return(u * sqrt( a * ( exp(- 2.0 / a * log(w)) - 1.0) / w)); 00086 } 00087 00088 void RandStudentT::shootArray( const int size, double* vect, 00089 double a ) 00090 { 00091 int i; 00092 00093 for (i=0; i<size; ++i) 00094 vect[i] = shoot(a); 00095 } 00096 00097 void RandStudentT::shootArray( HepRandomEngine* anEngine, 00098 const int size, double* vect, 00099 double a ) 00100 { 00101 int i; 00102 00103 for (i=0; i<size; ++i) 00104 vect[i] = shoot(anEngine,a); 00105 } 00106 00107 double RandStudentT::fire( double a ) { 00108 00109 double u,v,w; 00110 00111 do 00112 { 00113 u = 2.0 * localEngine->flat() - 1.0; 00114 v = 2.0 * localEngine->flat() - 1.0; 00115 } 00116 while ((w = u * u + v * v) > 1.0); 00117 00118 return(u * sqrt( a * ( exp(- 2.0 / a * log(w)) - 1.0) / w)); 00119 } 00120 00121 void RandStudentT::fireArray( const int size, double* vect) 00122 { 00123 int i; 00124 00125 for (i=0; i<size; ++i) 00126 vect[i] = fire(defaultA); 00127 } 00128 00129 void RandStudentT::fireArray( const int size, double* vect, 00130 double a ) 00131 { 00132 int i; 00133 00134 for (i=0; i<size; ++i) 00135 vect[i] = fire(a); 00136 } 00137 00138 double RandStudentT::shoot( HepRandomEngine *anEngine, double a ) { 00139 00140 double u,v,w; 00141 00142 do 00143 { 00144 u = 2.0 * anEngine->flat() - 1.0; 00145 v = 2.0 * anEngine->flat() - 1.0; 00146 } 00147 while ((w = u * u + v * v) > 1.0); 00148 00149 return(u * sqrt( a * ( exp(- 2.0 / a * log(w)) - 1.0) / w)); 00150 } 00151 00152 std::ostream & RandStudentT::put ( std::ostream & os ) const { 00153 int pr=os.precision(20); 00154 std::vector<unsigned long> t(2); 00155 os << " " << name() << "\n"; 00156 os << "Uvec" << "\n"; 00157 t = DoubConv::dto2longs(defaultA); 00158 os << defaultA << " " << t[0] << " " << t[1] << "\n"; 00159 os.precision(pr); 00160 return os; 00161 #ifdef REMOVED 00162 int pr=os.precision(20); 00163 os << " " << name() << "\n"; 00164 os << defaultA << "\n"; 00165 os.precision(pr); 00166 return os; 00167 #endif 00168 } 00169 00170 std::istream & RandStudentT::get ( std::istream & is ) { 00171 std::string inName; 00172 is >> inName; 00173 if (inName != name()) { 00174 is.clear(std::ios::badbit | is.rdstate()); 00175 std::cerr << "Mismatch when expecting to read state of a " 00176 << name() << " distribution\n" 00177 << "Name found was " << inName 00178 << "\nistream is left in the badbit state\n"; 00179 return is; 00180 } 00181 if (possibleKeywordInput(is, "Uvec", defaultA)) { 00182 std::vector<unsigned long> t(2); 00183 is >> defaultA >> t[0] >> t[1]; defaultA = DoubConv::longs2double(t); 00184 return is; 00185 } 00186 // is >> defaultA encompassed by possibleKeywordInput 00187 return is; 00188 } 00189 00190 } // namespace CLHEP