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