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