Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/externals/clhep/src/RandChiSquare.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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