CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: RandChiSquare.cc,v 1.6 2010/06/16 17:24:53 garren Exp $ 00002 // -*- C++ -*- 00003 // 00004 // ----------------------------------------------------------------------- 00005 // HEP Random 00006 // --- RandChiSquare --- 00007 // class implementation file 00008 // ----------------------------------------------------------------------- 00009 00010 // ======================================================================= 00011 // John Marraffino - Created: 12th May 1998 00012 // M Fischler - put and get to/from streams 12/10/04 00013 // M Fischler - put/get to/from streams uses pairs of ulongs when 00014 // + storing doubles avoid problems with precision 00015 // 4/14/05 00016 // ======================================================================= 00017 00018 #include "CLHEP/Random/defs.h" 00019 #include "CLHEP/Random/RandChiSquare.h" 00020 #include "CLHEP/Random/DoubConv.hh" 00021 #include <cmath> // for std::log() 00022 00023 namespace CLHEP { 00024 00025 std::string RandChiSquare::name() const {return "RandChiSquare";} 00026 HepRandomEngine & RandChiSquare::engine() {return *localEngine;} 00027 00028 RandChiSquare::~RandChiSquare() { 00029 } 00030 00031 double RandChiSquare::shoot( HepRandomEngine *anEngine, double a ) { 00032 return genChiSquare( anEngine, a ); 00033 } 00034 00035 double RandChiSquare::shoot( double a ) { 00036 HepRandomEngine *anEngine = HepRandom::getTheEngine(); 00037 return genChiSquare( anEngine, a ); 00038 } 00039 00040 double RandChiSquare::fire( double a ) { 00041 return genChiSquare( localEngine.get(), a ); 00042 } 00043 00044 void RandChiSquare::shootArray( const int size, double* vect, 00045 double a ) { 00046 for( double* v = vect; v != vect+size; ++v ) 00047 *v = shoot(a); 00048 } 00049 00050 void RandChiSquare::shootArray( HepRandomEngine* anEngine, 00051 const int size, double* vect, 00052 double a ) 00053 { 00054 for( double* v = vect; v != vect+size; ++v ) 00055 *v = shoot(anEngine,a); 00056 } 00057 00058 void RandChiSquare::fireArray( const int size, double* vect) { 00059 for( double* v = vect; v != vect+size; ++v ) 00060 *v = fire(defaultA); 00061 } 00062 00063 void RandChiSquare::fireArray( const int size, double* vect, 00064 double a ) { 00065 for( double* v = vect; v != vect+size; ++v ) 00066 *v = fire(a); 00067 } 00068 00069 double RandChiSquare::genChiSquare( HepRandomEngine *anEngine, 00070 double a ) { 00071 /****************************************************************** 00072 * * 00073 * Chi Distribution - Ratio of Uniforms with shift * 00074 * * 00075 ****************************************************************** 00076 * * 00077 * FUNCTION : - chru samples a random number from the Chi * 00078 * distribution with parameter a > 1. * 00079 * REFERENCE : - J.F. Monahan (1987): An algorithm for * 00080 * generating chi random variables, ACM Trans. * 00081 * Math. Software 13, 168-172. * 00082 * SUBPROGRAM : - anEngine ... pointer to a (0,1)-Uniform * 00083 * engine * 00084 * * 00085 * Implemented by R. Kremer, 1990 * 00086 ******************************************************************/ 00087 00088 static double a_in = -1.0,b,vm,vp,vd; 00089 double u,v,z,zz,r; 00090 00091 // Check for invalid input value 00092 00093 if( a < 1 ) return (-1.0); 00094 00095 if (a == 1) 00096 { 00097 for(;;) 00098 { 00099 u = anEngine->flat(); 00100 v = anEngine->flat() * 0.857763884960707; 00101 z = v / u; 00102 if (z < 0) continue; 00103 zz = z * z; 00104 r = 2.5 - zz; 00105 if (z < 0.0) r = r + zz * z / (3.0 * z); 00106 if (u < r * 0.3894003915) return(z*z); 00107 if (zz > (1.036961043 / u + 1.4)) continue; 00108 if (2 * std::log(u) < (- zz * 0.5 )) return(z*z); 00109 } 00110 } 00111 else 00112 { 00113 if (a != a_in) 00114 { 00115 b = std::sqrt(a - 1.0); 00116 vm = - 0.6065306597 * (1.0 - 0.25 / (b * b + 1.0)); 00117 vm = (-b > vm)? -b : vm; 00118 vp = 0.6065306597 * (0.7071067812 + b) / (0.5 + b); 00119 vd = vp - vm; 00120 a_in = a; 00121 } 00122 for(;;) 00123 { 00124 u = anEngine->flat(); 00125 v = anEngine->flat() * vd + vm; 00126 z = v / u; 00127 if (z < -b) continue; 00128 zz = z * z; 00129 r = 2.5 - zz; 00130 if (z < 0.0) r = r + zz * z / (3.0 * (z + b)); 00131 if (u < r * 0.3894003915) return((z + b)*(z + b)); 00132 if (zz > (1.036961043 / u + 1.4)) continue; 00133 if (2 * std::log(u) < (std::log(1.0 + z / b) * b * b - zz * 0.5 - z * b)) return((z + b)*(z + b)); 00134 } 00135 } 00136 } 00137 00138 std::ostream & RandChiSquare::put ( std::ostream & os ) const { 00139 int pr=os.precision(20); 00140 std::vector<unsigned long> t(2); 00141 os << " " << name() << "\n"; 00142 os << "Uvec" << "\n"; 00143 t = DoubConv::dto2longs(defaultA); 00144 os << defaultA << " " << t[0] << " " << t[1] << "\n"; 00145 os.precision(pr); 00146 return os; 00147 #ifdef REMOVED 00148 int pr=os.precision(20); 00149 os << " " << name() << "\n"; 00150 os << defaultA << "\n"; 00151 os.precision(pr); 00152 return os; 00153 #endif 00154 } 00155 00156 std::istream & RandChiSquare::get ( std::istream & is ) { 00157 std::string inName; 00158 is >> inName; 00159 if (inName != name()) { 00160 is.clear(std::ios::badbit | is.rdstate()); 00161 std::cerr << "Mismatch when expecting to read state of a " 00162 << name() << " distribution\n" 00163 << "Name found was " << inName 00164 << "\nistream is left in the badbit state\n"; 00165 return is; 00166 } 00167 if (possibleKeywordInput(is, "Uvec", defaultA)) { 00168 std::vector<unsigned long> t(2); 00169 is >> defaultA >> t[0] >> t[1]; defaultA = DoubConv::longs2double(t); 00170 return is; 00171 } 00172 // is >> defaultA encompassed by possibleKeywordInput 00173 return is; 00174 } 00175 00176 00177 00178 } // namespace CLHEP 00179