Geant4 Cross Reference |
>> 1 // $Id:$ 1 // -*- C++ -*- 2 // -*- C++ -*- 2 // 3 // 3 // ------------------------------------------- 4 // ----------------------------------------------------------------------- 4 // HEP Random 5 // HEP Random 5 // --- RandChiSquare 6 // --- RandChiSquare --- 6 // class implementation f 7 // class implementation file 7 // ------------------------------------------- 8 // ----------------------------------------------------------------------- 8 9 9 // =========================================== 10 // ======================================================================= 10 // John Marraffino - Created: 12th May 1998 11 // John Marraffino - Created: 12th May 1998 11 // M Fischler - put and get to/from stream 12 // M Fischler - put and get to/from streams 12/10/04 12 // M Fischler - put/get to/from streams 13 // M Fischler - put/get to/from streams uses pairs of ulongs when 13 // + storing doubles avoid problems with 14 // + storing doubles avoid problems with precision 14 // 4/14/05 15 // 4/14/05 15 // =========================================== 16 // ======================================================================= 16 17 17 #include "CLHEP/Random/RandChiSquare.h" 18 #include "CLHEP/Random/RandChiSquare.h" 18 #include "CLHEP/Random/DoubConv.h" 19 #include "CLHEP/Random/DoubConv.h" 19 #include "CLHEP/Utility/thread_local.h" 20 #include "CLHEP/Utility/thread_local.h" 20 #include <cmath> // for std::log() 21 #include <cmath> // for std::log() 21 #include <iostream> << 22 #include <string> << 23 #include <vector> << 24 22 25 namespace CLHEP { 23 namespace CLHEP { 26 24 27 std::string RandChiSquare::name() const {retur 25 std::string RandChiSquare::name() const {return "RandChiSquare";} 28 HepRandomEngine & RandChiSquare::engine() {ret 26 HepRandomEngine & RandChiSquare::engine() {return *localEngine;} 29 27 30 RandChiSquare::~RandChiSquare() { 28 RandChiSquare::~RandChiSquare() { 31 } 29 } 32 30 33 double RandChiSquare::shoot( HepRandomEngine * 31 double RandChiSquare::shoot( HepRandomEngine *anEngine, double a ) { 34 return genChiSquare( anEngine, a ); 32 return genChiSquare( anEngine, a ); 35 } 33 } 36 34 37 double RandChiSquare::shoot( double a ) { 35 double RandChiSquare::shoot( double a ) { 38 HepRandomEngine *anEngine = HepRandom::getTh 36 HepRandomEngine *anEngine = HepRandom::getTheEngine(); 39 return genChiSquare( anEngine, a ); 37 return genChiSquare( anEngine, a ); 40 } 38 } 41 39 42 double RandChiSquare::fire( double a ) { 40 double RandChiSquare::fire( double a ) { 43 return genChiSquare( localEngine.get(), a ); 41 return genChiSquare( localEngine.get(), a ); 44 } 42 } 45 43 46 void RandChiSquare::shootArray( const int size 44 void RandChiSquare::shootArray( const int size, double* vect, 47 double a ) { 45 double a ) { 48 for( double* v = vect; v != vect+size; ++v ) 46 for( double* v = vect; v != vect+size; ++v ) 49 *v = shoot(a); 47 *v = shoot(a); 50 } 48 } 51 49 52 void RandChiSquare::shootArray( HepRandomEngin 50 void RandChiSquare::shootArray( HepRandomEngine* anEngine, 53 const int size, do 51 const int size, double* vect, 54 double a ) 52 double a ) 55 { 53 { 56 for( double* v = vect; v != vect+size; ++v ) 54 for( double* v = vect; v != vect+size; ++v ) 57 *v = shoot(anEngine,a); 55 *v = shoot(anEngine,a); 58 } 56 } 59 57 60 void RandChiSquare::fireArray( const int size, 58 void RandChiSquare::fireArray( const int size, double* vect) { 61 for( double* v = vect; v != vect+size; ++v ) 59 for( double* v = vect; v != vect+size; ++v ) 62 *v = fire(defaultA); 60 *v = fire(defaultA); 63 } 61 } 64 62 65 void RandChiSquare::fireArray( const int size, 63 void RandChiSquare::fireArray( const int size, double* vect, 66 double a ) { 64 double a ) { 67 for( double* v = vect; v != vect+size; ++v ) 65 for( double* v = vect; v != vect+size; ++v ) 68 *v = fire(a); 66 *v = fire(a); 69 } 67 } 70 68 71 double RandChiSquare::genChiSquare( HepRandomE 69 double RandChiSquare::genChiSquare( HepRandomEngine *anEngine, 72 double 70 double a ) { 73 /********************************************* 71 /****************************************************************** 74 * 72 * * 75 * Chi Distribution - Ratio of Uniforms 73 * Chi Distribution - Ratio of Uniforms with shift * 76 * 74 * * 77 ********************************************* 75 ****************************************************************** 78 * 76 * * 79 * FUNCTION : - chru samples a random number 77 * FUNCTION : - chru samples a random number from the Chi * 80 * distribution with parameter 78 * distribution with parameter a > 1. * 81 * REFERENCE : - J.F. Monahan (1987): An algo 79 * REFERENCE : - J.F. Monahan (1987): An algorithm for * 82 * generating chi random variab 80 * generating chi random variables, ACM Trans. * 83 * Math. Software 13, 168-172. 81 * Math. Software 13, 168-172. * 84 * SUBPROGRAM : - anEngine ... pointer to a ( 82 * SUBPROGRAM : - anEngine ... pointer to a (0,1)-Uniform * 85 * engine 83 * engine * 86 * 84 * * 87 * Implemented by R. Kremer, 1990 85 * Implemented by R. Kremer, 1990 * 88 ********************************************* 86 ******************************************************************/ 89 87 90 static CLHEP_THREAD_LOCAL double a_in = -1.0, 88 static CLHEP_THREAD_LOCAL double a_in = -1.0,b,vm,vp,vd; 91 double u,v,z,zz,r; 89 double u,v,z,zz,r; 92 90 93 // Check for invalid input value 91 // Check for invalid input value 94 92 95 if( a < 1 ) return (-1.0); 93 if( a < 1 ) return (-1.0); 96 94 97 if (a == 1) 95 if (a == 1) 98 { 96 { 99 for(;;) 97 for(;;) 100 { 98 { 101 u = anEngine->flat(); 99 u = anEngine->flat(); 102 v = anEngine->flat() * 0.857763884960707; 100 v = anEngine->flat() * 0.857763884960707; 103 z = v / u; 101 z = v / u; 104 if (z < 0) continue; 102 if (z < 0) continue; 105 zz = z * z; 103 zz = z * z; 106 r = 2.5 - zz; 104 r = 2.5 - zz; 107 if (z < 0.0) r = r + zz * z / (3.0 * z); 105 if (z < 0.0) r = r + zz * z / (3.0 * z); 108 if (u < r * 0.3894003915) return(z*z); 106 if (u < r * 0.3894003915) return(z*z); 109 if (zz > (1.036961043 / u + 1.4)) continu 107 if (zz > (1.036961043 / u + 1.4)) continue; 110 if (2 * std::log(u) < (- zz * 0.5 )) retu 108 if (2 * std::log(u) < (- zz * 0.5 )) return(z*z); 111 } 109 } 112 } 110 } 113 else 111 else 114 { 112 { 115 if (a != a_in) 113 if (a != a_in) 116 { 114 { 117 b = std::sqrt(a - 1.0); 115 b = std::sqrt(a - 1.0); 118 vm = - 0.6065306597 * (1.0 - 0.25 / (b * 116 vm = - 0.6065306597 * (1.0 - 0.25 / (b * b + 1.0)); 119 vm = (-b > vm)? -b : vm; 117 vm = (-b > vm)? -b : vm; 120 vp = 0.6065306597 * (0.7071067812 + b) / 118 vp = 0.6065306597 * (0.7071067812 + b) / (0.5 + b); 121 vd = vp - vm; 119 vd = vp - vm; 122 a_in = a; 120 a_in = a; 123 } 121 } 124 for(;;) 122 for(;;) 125 { 123 { 126 u = anEngine->flat(); 124 u = anEngine->flat(); 127 v = anEngine->flat() * vd + vm; 125 v = anEngine->flat() * vd + vm; 128 z = v / u; 126 z = v / u; 129 if (z < -b) continue; 127 if (z < -b) continue; 130 zz = z * z; 128 zz = z * z; 131 r = 2.5 - zz; 129 r = 2.5 - zz; 132 if (z < 0.0) r = r + zz * z / (3.0 * (z + 130 if (z < 0.0) r = r + zz * z / (3.0 * (z + b)); 133 if (u < r * 0.3894003915) return((z + b)* 131 if (u < r * 0.3894003915) return((z + b)*(z + b)); 134 if (zz > (1.036961043 / u + 1.4)) continu 132 if (zz > (1.036961043 / u + 1.4)) continue; 135 if (2 * std::log(u) < (std::log(1.0 + z / 133 if (2 * std::log(u) < (std::log(1.0 + z / b) * b * b - zz * 0.5 - z * b)) return((z + b)*(z + b)); 136 } 134 } 137 } 135 } 138 } 136 } 139 137 140 std::ostream & RandChiSquare::put ( std::ostre 138 std::ostream & RandChiSquare::put ( std::ostream & os ) const { 141 long pr=os.precision(20); << 139 int pr=os.precision(20); 142 std::vector<unsigned long> t(2); 140 std::vector<unsigned long> t(2); 143 os << " " << name() << "\n"; 141 os << " " << name() << "\n"; 144 os << "Uvec" << "\n"; 142 os << "Uvec" << "\n"; 145 t = DoubConv::dto2longs(defaultA); 143 t = DoubConv::dto2longs(defaultA); 146 os << defaultA << " " << t[0] << " " << t[1] 144 os << defaultA << " " << t[0] << " " << t[1] << "\n"; 147 os.precision(pr); 145 os.precision(pr); 148 return os; 146 return os; 149 } 147 } 150 148 151 std::istream & RandChiSquare::get ( std::istre 149 std::istream & RandChiSquare::get ( std::istream & is ) { 152 std::string inName; 150 std::string inName; 153 is >> inName; 151 is >> inName; 154 if (inName != name()) { 152 if (inName != name()) { 155 is.clear(std::ios::badbit | is.rdstate()); 153 is.clear(std::ios::badbit | is.rdstate()); 156 std::cerr << "Mismatch when expecting to r 154 std::cerr << "Mismatch when expecting to read state of a " 157 << name() << " distribution\n" 155 << name() << " distribution\n" 158 << "Name found was " << inName 156 << "Name found was " << inName 159 << "\nistream is left in the badbit st 157 << "\nistream is left in the badbit state\n"; 160 return is; 158 return is; 161 } 159 } 162 if (possibleKeywordInput(is, "Uvec", default 160 if (possibleKeywordInput(is, "Uvec", defaultA)) { 163 std::vector<unsigned long> t(2); 161 std::vector<unsigned long> t(2); 164 is >> defaultA >> t[0] >> t[1]; defaultA = 162 is >> defaultA >> t[0] >> t[1]; defaultA = DoubConv::longs2double(t); 165 return is; 163 return is; 166 } 164 } 167 // is >> defaultA encompassed by possibleKey 165 // is >> defaultA encompassed by possibleKeywordInput 168 return is; 166 return is; 169 } 167 } 170 168 171 169 172 170 173 } // namespace CLHEP 171 } // namespace CLHEP 174 172 175 173