CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: RandChiSquare.cc,v 1.4.4.2 2005/04/15 16:32: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 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 if ( deleteEngine ) delete localEngine; 00030 } 00031 00032 RandChiSquare::RandChiSquare(const RandChiSquare& right) 00033 : defaultA(right.defaultA) 00034 {;} 00035 00036 double RandChiSquare::shoot( HepRandomEngine *anEngine, double a ) { 00037 return genChiSquare( anEngine, a ); 00038 } 00039 00040 double RandChiSquare::shoot( double a ) { 00041 HepRandomEngine *anEngine = HepRandom::getTheEngine(); 00042 return genChiSquare( anEngine, a ); 00043 } 00044 00045 double RandChiSquare::fire( double a ) { 00046 return genChiSquare( localEngine, a ); 00047 } 00048 00049 void RandChiSquare::shootArray( const int size, double* vect, 00050 double a ) { 00051 int i; 00052 00053 for (i=0; i<size; ++i) 00054 vect[i] = shoot(a); 00055 } 00056 00057 void RandChiSquare::shootArray( HepRandomEngine* anEngine, 00058 const int size, double* vect, 00059 double a ) 00060 { 00061 int i; 00062 00063 for (i=0; i<size; ++i) 00064 vect[i] = shoot(anEngine,a); 00065 } 00066 00067 void RandChiSquare::fireArray( const int size, double* vect) { 00068 int i; 00069 00070 for (i=0; i<size; ++i) 00071 vect[i] = fire(defaultA); 00072 } 00073 00074 void RandChiSquare::fireArray( const int size, double* vect, 00075 double a ) { 00076 int i; 00077 00078 for (i=0; i<size; ++i) 00079 vect[i] = fire(a); 00080 } 00081 00082 double RandChiSquare::genChiSquare( HepRandomEngine *anEngine, 00083 double a ) { 00084 /****************************************************************** 00085 * * 00086 * Chi Distribution - Ratio of Uniforms with shift * 00087 * * 00088 ****************************************************************** 00089 * * 00090 * FUNCTION : - chru samples a random number from the Chi * 00091 * distribution with parameter a > 1. * 00092 * REFERENCE : - J.F. Monahan (1987): An algorithm for * 00093 * generating chi random variables, ACM Trans. * 00094 * Math. Software 13, 168-172. * 00095 * SUBPROGRAM : - anEngine ... pointer to a (0,1)-Uniform * 00096 * engine * 00097 * * 00098 * Implemented by R. Kremer, 1990 * 00099 ******************************************************************/ 00100 00101 static double a_in = -1.0,b,vm,vp,vd; 00102 double u,v,z,zz,r; 00103 00104 // Check for invalid input value 00105 00106 if( a < 1 ) return (-1.0); 00107 00108 if (a == 1) 00109 { 00110 for(;;) 00111 { 00112 u = anEngine->flat(); 00113 v = anEngine->flat() * 0.857763884960707; 00114 z = v / u; 00115 if (z < 0) continue; 00116 zz = z * z; 00117 r = 2.5 - zz; 00118 if (z < 0.0) r = r + zz * z / (3.0 * z); 00119 if (u < r * 0.3894003915) return(z*z); 00120 if (zz > (1.036961043 / u + 1.4)) continue; 00121 if (2 * log(u) < (- zz * 0.5 )) return(z*z); 00122 } 00123 } 00124 else 00125 { 00126 if (a != a_in) 00127 { 00128 b = sqrt(a - 1.0); 00129 vm = - 0.6065306597 * (1.0 - 0.25 / (b * b + 1.0)); 00130 vm = (-b > vm)? -b : vm; 00131 vp = 0.6065306597 * (0.7071067812 + b) / (0.5 + b); 00132 vd = vp - vm; 00133 a_in = a; 00134 } 00135 for(;;) 00136 { 00137 u = anEngine->flat(); 00138 v = anEngine->flat() * vd + vm; 00139 z = v / u; 00140 if (z < -b) continue; 00141 zz = z * z; 00142 r = 2.5 - zz; 00143 if (z < 0.0) r = r + zz * z / (3.0 * (z + b)); 00144 if (u < r * 0.3894003915) return((z + b)*(z + b)); 00145 if (zz > (1.036961043 / u + 1.4)) continue; 00146 if (2 * log(u) < (log(1.0 + z / b) * b * b - zz * 0.5 - z * b)) return((z + b)*(z + b)); 00147 } 00148 } 00149 } 00150 00151 std::ostream & RandChiSquare::put ( std::ostream & os ) const { 00152 int pr=os.precision(20); 00153 std::vector<unsigned long> t(2); 00154 os << " " << name() << "\n"; 00155 os << "Uvec" << "\n"; 00156 t = DoubConv::dto2longs(defaultA); 00157 os << defaultA << " " << t[0] << " " << t[1] << "\n"; 00158 os.precision(pr); 00159 return os; 00160 #ifdef REMOVED 00161 int pr=os.precision(20); 00162 os << " " << name() << "\n"; 00163 os << defaultA << "\n"; 00164 os.precision(pr); 00165 return os; 00166 #endif 00167 } 00168 00169 std::istream & RandChiSquare::get ( std::istream & is ) { 00170 std::string inName; 00171 is >> inName; 00172 if (inName != name()) { 00173 is.clear(std::ios::badbit | is.rdstate()); 00174 std::cerr << "Mismatch when expecting to read state of a " 00175 << name() << " distribution\n" 00176 << "Name found was " << inName 00177 << "\nistream is left in the badbit state\n"; 00178 return is; 00179 } 00180 if (possibleKeywordInput(is, "Uvec", defaultA)) { 00181 std::vector<unsigned long> t(2); 00182 is >> defaultA >> t[0] >> t[1]; defaultA = DoubConv::longs2double(t); 00183 return is; 00184 } 00185 // is >> defaultA encompassed by possibleKeywordInput 00186 return is; 00187 } 00188 00189 00190 00191 } // namespace CLHEP 00192