Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/externals/clhep/src/RandStudentT.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 //                          --- RandStudentT ---
  6 //                      class implementation file
  7 // -----------------------------------------------------------------------
  8 
  9 // =======================================================================
 10 // John Marraffino - Created: 12th May 1998
 11 // G.Cosmo         - Fixed minor bug on inline definition for shoot()
 12 //                   methods : 20th Aug 1998
 13 // M Fischler      - put and get to/from streams 12/13/04
 14 // M Fischler       - put/get to/from streams uses pairs of ulongs when
 15 //      + storing doubles avoid problems with precision 
 16 //      4/14/05
 17 // =======================================================================
 18 
 19 #include <float.h>
 20 #include "CLHEP/Random/RandStudentT.h"
 21 #include "CLHEP/Random/DoubConv.h"
 22 
 23 #include <cmath>  // for std::log() std::exp()
 24 #include <iostream>
 25 #include <string>
 26 #include <vector>
 27 
 28 namespace CLHEP {
 29 
 30 std::string RandStudentT::name() const {return "RandStudentT";}
 31 HepRandomEngine & RandStudentT::engine() {return *localEngine;}
 32 
 33 RandStudentT::~RandStudentT()
 34 {
 35 }
 36 
 37 double RandStudentT::operator()() {
 38   return fire( defaultA );
 39 }
 40 
 41 double RandStudentT::operator()( double a ) {
 42   return fire( a );
 43 }
 44 
 45 double RandStudentT::shoot( double a ) {
 46 /******************************************************************
 47  *                                                                *
 48  *           Student-t Distribution - Polar Method                *
 49  *                                                                *
 50  ******************************************************************
 51  *                                                                *
 52  * The polar method of Box/Muller for generating Normal variates  *
 53  * is adapted to the Student-t distribution. The two generated    *
 54  * variates are not independent and the expected no. of uniforms  *
 55  * per variate is 2.5464.                                         *
 56  *                                                                *
 57  ******************************************************************
 58  *                                                                *
 59  * FUNCTION :   - tpol  samples a random number from the          *
 60  *                Student-t distribution with parameter a > 0.    *
 61  * REFERENCE :  - R.W. Bailey (1994): Polar generation of random  *
 62  *                variates with the t-distribution, Mathematics   *
 63  *                of Computation 62, 779-781.                     *
 64  * SUBPROGRAM : -  ... (0,1)-Uniform generator                    *
 65  *                                                                *
 66  *                                                                *
 67  * Implemented by F. Niederl, 1994                                *
 68  ******************************************************************/
 69 
 70  double u,v,w;
 71 
 72 // Check for valid input value
 73 
 74  if( a < 0.0) return (DBL_MAX);
 75 
 76  do
 77  {
 78    u = 2.0 * HepRandom::getTheEngine()->flat() - 1.0;
 79    v = 2.0 * HepRandom::getTheEngine()->flat() - 1.0;
 80  }
 81  while ((w = u * u + v * v) > 1.0);
 82 
 83  return(u * std::sqrt( a * ( std::exp(- 2.0 / a * std::log(w)) - 1.0) / w));
 84 }
 85 
 86 void RandStudentT::shootArray( const int size, double* vect,
 87                             double a )
 88 {
 89   for( double* v = vect; v != vect + size; ++v )
 90     *v = shoot(a);
 91 }
 92 
 93 void RandStudentT::shootArray( HepRandomEngine* anEngine,
 94                             const int size, double* vect,
 95                             double a )
 96 {
 97   for( double* v = vect; v != vect + size; ++v )
 98     *v = shoot(anEngine,a);
 99 }
100 
101 double RandStudentT::fire( double a ) {
102 
103  double u,v,w;
104 
105  do
106  {
107    u = 2.0 * localEngine->flat() - 1.0;
108    v = 2.0 * localEngine->flat() - 1.0;
109  }
110  while ((w = u * u + v * v) > 1.0);
111 
112  return(u * std::sqrt( a * ( std::exp(- 2.0 / a * std::log(w)) - 1.0) / w));
113 }
114 
115 void RandStudentT::fireArray( const int size, double* vect)
116 {
117   for( double* v = vect; v != vect + size; ++v )
118     *v = fire(defaultA);
119 }
120 
121 void RandStudentT::fireArray( const int size, double* vect,
122                            double a )
123 {
124   for( double* v = vect; v != vect + size; ++v )
125     *v = fire(a);
126 }
127 
128 double RandStudentT::shoot( HepRandomEngine *anEngine, double a ) {
129 
130  double u,v,w;
131 
132  do
133  {
134    u = 2.0 * anEngine->flat() - 1.0;
135    v = 2.0 * anEngine->flat() - 1.0;
136  }
137  while ((w = u * u + v * v) > 1.0);
138 
139  return(u * std::sqrt( a * ( std::exp(- 2.0 / a * std::log(w)) - 1.0) / w));
140 }
141 
142 std::ostream & RandStudentT::put ( std::ostream & os ) const {
143   long pr=os.precision(20);
144   std::vector<unsigned long> t(2);
145   os << " " << name() << "\n";
146   os << "Uvec" << "\n";
147   t = DoubConv::dto2longs(defaultA);
148   os << defaultA << " " << t[0] << " " << t[1] << "\n";
149   os.precision(pr);
150   return os;
151 }
152 
153 std::istream & RandStudentT::get ( std::istream & is ) {
154   std::string inName;
155   is >> inName;
156   if (inName != name()) {
157     is.clear(std::ios::badbit | is.rdstate());
158     std::cerr << "Mismatch when expecting to read state of a "
159             << name() << " distribution\n"
160         << "Name found was " << inName
161         << "\nistream is left in the badbit state\n";
162     return is;
163   }
164   if (possibleKeywordInput(is, "Uvec", defaultA)) {
165     std::vector<unsigned long> t(2);
166     is >> defaultA >> t[0] >> t[1]; defaultA = DoubConv::longs2double(t); 
167     return is;
168   }
169   // is >> defaultA encompassed by possibleKeywordInput
170   return is;
171 }
172 
173 }  // namespace CLHEP
174