00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include "CLHEP/Random/RandChiSquare.h"
00019 #include "CLHEP/Random/DoubConv.h"
00020 #include <cmath>
00021
00022 namespace CLHEP {
00023
00024 std::string RandChiSquare::name() const {return "RandChiSquare";}
00025 HepRandomEngine & RandChiSquare::engine() {return *localEngine;}
00026
00027 RandChiSquare::~RandChiSquare() {
00028 }
00029
00030 double RandChiSquare::shoot( HepRandomEngine *anEngine, double a ) {
00031 return genChiSquare( anEngine, a );
00032 }
00033
00034 double RandChiSquare::shoot( double a ) {
00035 HepRandomEngine *anEngine = HepRandom::getTheEngine();
00036 return genChiSquare( anEngine, a );
00037 }
00038
00039 double RandChiSquare::fire( double a ) {
00040 return genChiSquare( localEngine.get(), a );
00041 }
00042
00043 void RandChiSquare::shootArray( const int size, double* vect,
00044 double a ) {
00045 for( double* v = vect; v != vect+size; ++v )
00046 *v = shoot(a);
00047 }
00048
00049 void RandChiSquare::shootArray( HepRandomEngine* anEngine,
00050 const int size, double* vect,
00051 double a )
00052 {
00053 for( double* v = vect; v != vect+size; ++v )
00054 *v = shoot(anEngine,a);
00055 }
00056
00057 void RandChiSquare::fireArray( const int size, double* vect) {
00058 for( double* v = vect; v != vect+size; ++v )
00059 *v = fire(defaultA);
00060 }
00061
00062 void RandChiSquare::fireArray( const int size, double* vect,
00063 double a ) {
00064 for( double* v = vect; v != vect+size; ++v )
00065 *v = fire(a);
00066 }
00067
00068 double RandChiSquare::genChiSquare( HepRandomEngine *anEngine,
00069 double a ) {
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087 static double a_in = -1.0,b,vm,vp,vd;
00088 double u,v,z,zz,r;
00089
00090
00091
00092 if( a < 1 ) return (-1.0);
00093
00094 if (a == 1)
00095 {
00096 for(;;)
00097 {
00098 u = anEngine->flat();
00099 v = anEngine->flat() * 0.857763884960707;
00100 z = v / u;
00101 if (z < 0) continue;
00102 zz = z * z;
00103 r = 2.5 - zz;
00104 if (z < 0.0) r = r + zz * z / (3.0 * z);
00105 if (u < r * 0.3894003915) return(z*z);
00106 if (zz > (1.036961043 / u + 1.4)) continue;
00107 if (2 * std::log(u) < (- zz * 0.5 )) return(z*z);
00108 }
00109 }
00110 else
00111 {
00112 if (a != a_in)
00113 {
00114 b = std::sqrt(a - 1.0);
00115 vm = - 0.6065306597 * (1.0 - 0.25 / (b * b + 1.0));
00116 vm = (-b > vm)? -b : vm;
00117 vp = 0.6065306597 * (0.7071067812 + b) / (0.5 + b);
00118 vd = vp - vm;
00119 a_in = a;
00120 }
00121 for(;;)
00122 {
00123 u = anEngine->flat();
00124 v = anEngine->flat() * vd + vm;
00125 z = v / u;
00126 if (z < -b) continue;
00127 zz = z * z;
00128 r = 2.5 - zz;
00129 if (z < 0.0) r = r + zz * z / (3.0 * (z + b));
00130 if (u < r * 0.3894003915) return((z + b)*(z + b));
00131 if (zz > (1.036961043 / u + 1.4)) continue;
00132 if (2 * std::log(u) < (std::log(1.0 + z / b) * b * b - zz * 0.5 - z * b)) return((z + b)*(z + b));
00133 }
00134 }
00135 }
00136
00137 std::ostream & RandChiSquare::put ( std::ostream & os ) const {
00138 int pr=os.precision(20);
00139 std::vector<unsigned long> t(2);
00140 os << " " << name() << "\n";
00141 os << "Uvec" << "\n";
00142 t = DoubConv::dto2longs(defaultA);
00143 os << defaultA << " " << t[0] << " " << t[1] << "\n";
00144 os.precision(pr);
00145 return os;
00146 }
00147
00148 std::istream & RandChiSquare::get ( std::istream & is ) {
00149 std::string inName;
00150 is >> inName;
00151 if (inName != name()) {
00152 is.clear(std::ios::badbit | is.rdstate());
00153 std::cerr << "Mismatch when expecting to read state of a "
00154 << name() << " distribution\n"
00155 << "Name found was " << inName
00156 << "\nistream is left in the badbit state\n";
00157 return is;
00158 }
00159 if (possibleKeywordInput(is, "Uvec", defaultA)) {
00160 std::vector<unsigned long> t(2);
00161 is >> defaultA >> t[0] >> t[1]; defaultA = DoubConv::longs2double(t);
00162 return is;
00163 }
00164
00165 return is;
00166 }
00167
00168
00169
00170 }
00171