CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: RandStudentT.cc,v 1.6 2010/06/16 17:24: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 std::log() std::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 { 00033 } 00034 00035 double RandStudentT::operator()() { 00036 return fire( defaultA ); 00037 } 00038 00039 double RandStudentT::operator()( double a ) { 00040 return fire( a ); 00041 } 00042 00043 double RandStudentT::shoot( double a ) { 00044 /****************************************************************** 00045 * * 00046 * Student-t Distribution - Polar Method * 00047 * * 00048 ****************************************************************** 00049 * * 00050 * The polar method of Box/Muller for generating Normal variates * 00051 * is adapted to the Student-t distribution. The two generated * 00052 * variates are not independent and the expected no. of uniforms * 00053 * per variate is 2.5464. * 00054 * * 00055 ****************************************************************** 00056 * * 00057 * FUNCTION : - tpol samples a random number from the * 00058 * Student-t distribution with parameter a > 0. * 00059 * REFERENCE : - R.W. Bailey (1994): Polar generation of random * 00060 * variates with the t-distribution, Mathematics * 00061 * of Computation 62, 779-781. * 00062 * SUBPROGRAM : - ... (0,1)-Uniform generator * 00063 * * 00064 * * 00065 * Implemented by F. Niederl, 1994 * 00066 ******************************************************************/ 00067 00068 double u,v,w; 00069 00070 // Check for valid input value 00071 00072 if( a < 0.0) return (DBL_MAX); 00073 00074 do 00075 { 00076 u = 2.0 * HepRandom::getTheEngine()->flat() - 1.0; 00077 v = 2.0 * HepRandom::getTheEngine()->flat() - 1.0; 00078 } 00079 while ((w = u * u + v * v) > 1.0); 00080 00081 return(u * std::sqrt( a * ( std::exp(- 2.0 / a * std::log(w)) - 1.0) / w)); 00082 } 00083 00084 void RandStudentT::shootArray( const int size, double* vect, 00085 double a ) 00086 { 00087 for( double* v = vect; v != vect + size; ++v ) 00088 *v = shoot(a); 00089 } 00090 00091 void RandStudentT::shootArray( HepRandomEngine* anEngine, 00092 const int size, double* vect, 00093 double a ) 00094 { 00095 for( double* v = vect; v != vect + size; ++v ) 00096 *v = shoot(anEngine,a); 00097 } 00098 00099 double RandStudentT::fire( double a ) { 00100 00101 double u,v,w; 00102 00103 do 00104 { 00105 u = 2.0 * localEngine->flat() - 1.0; 00106 v = 2.0 * localEngine->flat() - 1.0; 00107 } 00108 while ((w = u * u + v * v) > 1.0); 00109 00110 return(u * std::sqrt( a * ( std::exp(- 2.0 / a * std::log(w)) - 1.0) / w)); 00111 } 00112 00113 void RandStudentT::fireArray( const int size, double* vect) 00114 { 00115 for( double* v = vect; v != vect + size; ++v ) 00116 *v = fire(defaultA); 00117 } 00118 00119 void RandStudentT::fireArray( const int size, double* vect, 00120 double a ) 00121 { 00122 for( double* v = vect; v != vect + size; ++v ) 00123 *v = fire(a); 00124 } 00125 00126 double RandStudentT::shoot( HepRandomEngine *anEngine, double a ) { 00127 00128 double u,v,w; 00129 00130 do 00131 { 00132 u = 2.0 * anEngine->flat() - 1.0; 00133 v = 2.0 * anEngine->flat() - 1.0; 00134 } 00135 while ((w = u * u + v * v) > 1.0); 00136 00137 return(u * std::sqrt( a * ( std::exp(- 2.0 / a * std::log(w)) - 1.0) / w)); 00138 } 00139 00140 std::ostream & RandStudentT::put ( std::ostream & os ) const { 00141 int pr=os.precision(20); 00142 std::vector<unsigned long> t(2); 00143 os << " " << name() << "\n"; 00144 os << "Uvec" << "\n"; 00145 t = DoubConv::dto2longs(defaultA); 00146 os << defaultA << " " << t[0] << " " << t[1] << "\n"; 00147 os.precision(pr); 00148 return os; 00149 #ifdef REMOVED 00150 int pr=os.precision(20); 00151 os << " " << name() << "\n"; 00152 os << defaultA << "\n"; 00153 os.precision(pr); 00154 return os; 00155 #endif 00156 } 00157 00158 std::istream & RandStudentT::get ( std::istream & is ) { 00159 std::string inName; 00160 is >> inName; 00161 if (inName != name()) { 00162 is.clear(std::ios::badbit | is.rdstate()); 00163 std::cerr << "Mismatch when expecting to read state of a " 00164 << name() << " distribution\n" 00165 << "Name found was " << inName 00166 << "\nistream is left in the badbit state\n"; 00167 return is; 00168 } 00169 if (possibleKeywordInput(is, "Uvec", defaultA)) { 00170 std::vector<unsigned long> t(2); 00171 is >> defaultA >> t[0] >> t[1]; defaultA = DoubConv::longs2double(t); 00172 return is; 00173 } 00174 // is >> defaultA encompassed by possibleKeywordInput 00175 return is; 00176 } 00177 00178 } // namespace CLHEP