| Geant4 Cross Reference |
1 // -*- C++ -*- 1 // -*- C++ -*-
2 // 2 //
3 // ------------------------------------------- 3 // -----------------------------------------------------------------------
4 // HEP Random 4 // HEP Random
5 // --- RandGeneral -- 5 // --- RandGeneral ---
6 // class implementation f 6 // class implementation file
7 // ------------------------------------------- 7 // -----------------------------------------------------------------------
8 8
9 // =========================================== 9 // =======================================================================
10 // S.Magni & G.Pieri - Created: 5th September 10 // S.Magni & G.Pieri - Created: 5th September 1995
11 // G.Cosmo - Added constructor using 11 // G.Cosmo - Added constructor using default engine from the
12 // static generator. Simpl 12 // static generator. Simplified shoot() and
13 // shootArray() (not neede 13 // shootArray() (not needed in principle!): 20th Aug 1998
14 // M.G.Pia & G.Cosmo - Fixed bug in computatio 14 // M.G.Pia & G.Cosmo - Fixed bug in computation of theIntegralPdf in
15 // two constructors: 5th J 15 // two constructors: 5th Jan 1999
16 // S.Magni & G.Pieri - Added linear interpolat 16 // S.Magni & G.Pieri - Added linear interpolation: 24th Mar 1999
17 // M. Fischler - General cleanup: 14th M 17 // M. Fischler - General cleanup: 14th May 1999
18 // + Eliminated constructor code replicat 18 // + Eliminated constructor code replication by factoring
19 // common code into prepareTable. 19 // common code into prepareTable.
20 // + Eliminated fire/shoot code replicati 20 // + Eliminated fire/shoot code replication by factoring
21 // out common code into mapRandom. 21 // out common code into mapRandom.
22 // + A couple of methods are moved inline 22 // + A couple of methods are moved inline to avoid a
23 // speed cost for factoring out mapRand 23 // speed cost for factoring out mapRandom: fire()
24 // and shoot(anEngine). 24 // and shoot(anEngine).
25 // + Inserted checks for negative weight 25 // + Inserted checks for negative weight and zero total
26 // weight in the bins. 26 // weight in the bins.
27 // + Modified the binary search loop to a 27 // + Modified the binary search loop to avoid incorrect
28 // behavior when rand is example on a b 28 // behavior when rand is example on a boundary.
29 // + Moved the check of InterpolationType 29 // + Moved the check of InterpolationType up into
30 // the constructor. A type other than 30 // the constructor. A type other than 0 or 1
31 // will give the interpolated distribut 31 // will give the interpolated distribution (instead of
32 // a distribution that always returns 0 32 // a distribution that always returns 0).
33 // + Modified the computation of the retu 33 // + Modified the computation of the returned value
34 // to use algeraic simplification to im 34 // to use algeraic simplification to improve speed.
35 // Eliminated two of the three division 35 // Eliminated two of the three divisionns, made
36 // use of the fact that nabove-nbelow i 36 // use of the fact that nabove-nbelow is always 1, etc.
37 // + Inserted a check for rand hitting th 37 // + Inserted a check for rand hitting the boundary of a
38 // zero-width bin, to avoid dividing 0/ 38 // zero-width bin, to avoid dividing 0/0.
39 // M. Fischler - Minor correction in as 39 // M. Fischler - Minor correction in assert 31 July 2001
40 // + changed from assert (above = below+1 40 // + changed from assert (above = below+1) to ==
41 // M Fischler - put and get to/from st 41 // M Fischler - put and get to/from streams 12/15/04
42 // + Modifications to use a vector as the 42 // + Modifications to use a vector as theIntegraPdf
43 // M Fischler - put/get to/from streams 43 // M Fischler - put/get to/from streams uses pairs of ulongs when
44 // + storing doubles avoid problems with 44 // + storing doubles avoid problems with precision
45 // 4/14/05 45 // 4/14/05
46 // 46 //
47 // =========================================== 47 // =======================================================================
48 48
49 #include "CLHEP/Random/RandGeneral.h" 49 #include "CLHEP/Random/RandGeneral.h"
50 #include "CLHEP/Random/DoubConv.h" 50 #include "CLHEP/Random/DoubConv.h"
51 #include <cassert> 51 #include <cassert>
52 #include <iostream> 52 #include <iostream>
53 #include <string> 53 #include <string>
54 #include <vector> 54 #include <vector>
55 55
56 namespace CLHEP { 56 namespace CLHEP {
57 57
58 std::string RandGeneral::name() const {return 58 std::string RandGeneral::name() const {return "RandGeneral";}
59 HepRandomEngine & RandGeneral::engine() {retur 59 HepRandomEngine & RandGeneral::engine() {return *localEngine;}
60 60
61 61
62 ////////////////// 62 //////////////////
63 // Constructors 63 // Constructors
64 ////////////////// 64 //////////////////
65 65
66 RandGeneral::RandGeneral( const double* aProbF 66 RandGeneral::RandGeneral( const double* aProbFunc,
67 int theProbSize, 67 int theProbSize,
68 int IntType ) 68 int IntType )
69 : HepRandom(), 69 : HepRandom(),
70 localEngine(HepRandom::getTheEngine(), do_ 70 localEngine(HepRandom::getTheEngine(), do_nothing_deleter()),
71 nBins(theProbSize), 71 nBins(theProbSize),
72 InterpolationType(IntType) 72 InterpolationType(IntType)
73 { 73 {
74 prepareTable(aProbFunc); 74 prepareTable(aProbFunc);
75 } 75 }
76 76
77 RandGeneral::RandGeneral(HepRandomEngine& anEn 77 RandGeneral::RandGeneral(HepRandomEngine& anEngine,
78 const double* aProbFu 78 const double* aProbFunc,
79 int theProbSize, 79 int theProbSize,
80 int IntType ) 80 int IntType )
81 : HepRandom(), 81 : HepRandom(),
82 localEngine(&anEngine, do_nothing_deleter()) 82 localEngine(&anEngine, do_nothing_deleter()),
83 nBins(theProbSize), 83 nBins(theProbSize),
84 InterpolationType(IntType) 84 InterpolationType(IntType)
85 { 85 {
86 prepareTable(aProbFunc); 86 prepareTable(aProbFunc);
87 } 87 }
88 88
89 RandGeneral::RandGeneral(HepRandomEngine* anEn 89 RandGeneral::RandGeneral(HepRandomEngine* anEngine,
90 const double* aProbFu 90 const double* aProbFunc,
91 int theProbSize, 91 int theProbSize,
92 int IntType ) 92 int IntType )
93 : HepRandom(), 93 : HepRandom(),
94 localEngine(anEngine), 94 localEngine(anEngine),
95 nBins(theProbSize), 95 nBins(theProbSize),
96 InterpolationType(IntType) 96 InterpolationType(IntType)
97 { 97 {
98 prepareTable(aProbFunc); 98 prepareTable(aProbFunc);
99 } 99 }
100 100
101 void RandGeneral::prepareTable(const double* a 101 void RandGeneral::prepareTable(const double* aProbFunc) {
102 // 102 //
103 // Private method called only by constructors. 103 // Private method called only by constructors. Prepares theIntegralPdf.
104 // 104 //
105 if (nBins < 1) { 105 if (nBins < 1) {
106 std::cerr << 106 std::cerr <<
107 "RandGeneral constructed with no bins - will 107 "RandGeneral constructed with no bins - will use flat distribution\n";
108 useFlatDistribution(); 108 useFlatDistribution();
109 return; 109 return;
110 } 110 }
111 111
112 theIntegralPdf.resize(nBins+1); 112 theIntegralPdf.resize(nBins+1);
113 theIntegralPdf[0] = 0; 113 theIntegralPdf[0] = 0;
114 int ptn; 114 int ptn;
115 double weight; 115 double weight;
116 116
117 for ( ptn = 0; ptn<nBins; ++ptn ) { 117 for ( ptn = 0; ptn<nBins; ++ptn ) {
118 weight = aProbFunc[ptn]; 118 weight = aProbFunc[ptn];
119 if ( weight < 0 ) { 119 if ( weight < 0 ) {
120 // We can't stomach negative bin contents, 120 // We can't stomach negative bin contents, they invalidate the
121 // search algorithm when the distribution 121 // search algorithm when the distribution is fired.
122 std::cerr << 122 std::cerr <<
123 "RandGeneral constructed with negative-weigh 123 "RandGeneral constructed with negative-weight bin " << ptn <<
124 " = " << weight << " \n -- will substitute 124 " = " << weight << " \n -- will substitute 0 weight \n";
125 weight = 0; 125 weight = 0;
126 } 126 }
127 // std::cout << ptn << " " << weight << " 127 // std::cout << ptn << " " << weight << " " << theIntegralPdf[ptn] << "\n";
128 theIntegralPdf[ptn+1] = theIntegralPdf[ptn 128 theIntegralPdf[ptn+1] = theIntegralPdf[ptn] + weight;
129 } 129 }
130 130
131 if ( theIntegralPdf[nBins] <= 0 ) { 131 if ( theIntegralPdf[nBins] <= 0 ) {
132 std::cerr << 132 std::cerr <<
133 "RandGeneral constructed nothing in bins 133 "RandGeneral constructed nothing in bins - will use flat distribution\n";
134 useFlatDistribution(); 134 useFlatDistribution();
135 return; 135 return;
136 } 136 }
137 137
138 for ( ptn = 0; ptn < nBins+1; ++ptn ) { 138 for ( ptn = 0; ptn < nBins+1; ++ptn ) {
139 theIntegralPdf[ptn] /= theIntegralPdf[nBin 139 theIntegralPdf[ptn] /= theIntegralPdf[nBins];
140 // std::cout << ptn << " " << theIntegral 140 // std::cout << ptn << " " << theIntegralPdf[ptn] << "\n";
141 } 141 }
142 142
143 // And another useful variable is ... 143 // And another useful variable is ...
144 oneOverNbins = 1.0 / nBins; 144 oneOverNbins = 1.0 / nBins;
145 145
146 // One last chore: 146 // One last chore:
147 147
148 if ( (InterpolationType != 0) && (Interpolat 148 if ( (InterpolationType != 0) && (InterpolationType != 1) ) {
149 std::cerr << 149 std::cerr <<
150 "RandGeneral does not recognize IntType 150 "RandGeneral does not recognize IntType " << InterpolationType
151 << "\n Will use type 0 (continuous linea 151 << "\n Will use type 0 (continuous linear interpolation \n";
152 InterpolationType = 0; 152 InterpolationType = 0;
153 } 153 }
154 154
155 } // prepareTable() 155 } // prepareTable()
156 156
157 void RandGeneral::useFlatDistribution() { 157 void RandGeneral::useFlatDistribution() {
158 // 158 //
159 // Private method called only by prepareTables 159 // Private method called only by prepareTables in case of user error.
160 // 160 //
161 nBins = 1; 161 nBins = 1;
162 theIntegralPdf.resize(2); 162 theIntegralPdf.resize(2);
163 theIntegralPdf[0] = 0; 163 theIntegralPdf[0] = 0;
164 theIntegralPdf[1] = 1; 164 theIntegralPdf[1] = 1;
165 oneOverNbins = 1.0; 165 oneOverNbins = 1.0;
166 return; 166 return;
167 167
168 } // UseFlatDistribution() 168 } // UseFlatDistribution()
169 169
170 ////////////////// 170 //////////////////
171 // Destructor 171 // Destructor
172 ////////////////// 172 //////////////////
173 173
174 RandGeneral::~RandGeneral() { 174 RandGeneral::~RandGeneral() {
175 } 175 }
176 176
177 177
178 /////////////////// 178 ///////////////////
179 // mapRandom(rand) 179 // mapRandom(rand)
180 /////////////////// 180 ///////////////////
181 181
182 double RandGeneral::mapRandom(double rand) con 182 double RandGeneral::mapRandom(double rand) const {
183 // 183 //
184 // Private method to take the random (however 184 // Private method to take the random (however it is created) and map it
185 // according to the distribution. 185 // according to the distribution.
186 // 186 //
187 187
188 int nbelow = 0; // largest k such that I[k 188 int nbelow = 0; // largest k such that I[k] is known to be <= rand
189 int nabove = nBins; // largest k such th 189 int nabove = nBins; // largest k such that I[k] is known to be > rand
190 int middle; 190 int middle;
191 191
192 while (nabove > nbelow+1) { 192 while (nabove > nbelow+1) {
193 middle = (nabove + nbelow+1)>>1; 193 middle = (nabove + nbelow+1)>>1;
194 if (rand >= theIntegralPdf[middle]) { 194 if (rand >= theIntegralPdf[middle]) {
195 nbelow = middle; 195 nbelow = middle;
196 } else { 196 } else {
197 nabove = middle; 197 nabove = middle;
198 } 198 }
199 } // after this loop, nabove is always nbelo 199 } // after this loop, nabove is always nbelow+1 and they straddle rad:
200 assert ( nabove == nbelow+1 ); 200 assert ( nabove == nbelow+1 );
201 assert ( theIntegralPdf[nbelow] <= rand ); 201 assert ( theIntegralPdf[nbelow] <= rand );
202 assert ( theIntegralPdf[nabove] >= rand ); 202 assert ( theIntegralPdf[nabove] >= rand );
203 // If a defective engine produces rand=1, 203 // If a defective engine produces rand=1, that will
204 // still give sensible results so we relax 204 // still give sensible results so we relax the > rand assertion
205 205
206 if ( InterpolationType == 1 ) { 206 if ( InterpolationType == 1 ) {
207 207
208 return nbelow * oneOverNbins; 208 return nbelow * oneOverNbins;
209 209
210 } else { 210 } else {
211 211
212 double binMeasure = theIntegralPdf[nabove] 212 double binMeasure = theIntegralPdf[nabove] - theIntegralPdf[nbelow];
213 // binMeasure is always aProbFunc[nbelow], 213 // binMeasure is always aProbFunc[nbelow],
214 // but we don't have aProbFunc any more so 214 // but we don't have aProbFunc any more so we subtract.
215 215
216 if ( binMeasure == 0 ) { 216 if ( binMeasure == 0 ) {
217 // rand lies right in a bin of measure 0. S 217 // rand lies right in a bin of measure 0. Simply return the center
218 // of the range of that bin. (Any value bet 218 // of the range of that bin. (Any value between k/N and (k+1)/N is
219 // equally good, in this rare case.) 219 // equally good, in this rare case.)
220 return (nbelow + .5) * oneOverNbins; 220 return (nbelow + .5) * oneOverNbins;
221 } 221 }
222 222
223 double binFraction = (rand - theIntegralPd 223 double binFraction = (rand - theIntegralPdf[nbelow]) / binMeasure;
224 224
225 return (nbelow + binFraction) * oneOverNbi 225 return (nbelow + binFraction) * oneOverNbins;
226 } 226 }
227 227
228 } // mapRandom(rand) 228 } // mapRandom(rand)
229 229
230 void RandGeneral::shootArray( HepRandomEngine* 230 void RandGeneral::shootArray( HepRandomEngine* anEngine,
231 const int size, do 231 const int size, double* vect )
232 { 232 {
233 int i; 233 int i;
234 234
235 for (i=0; i<size; ++i) { 235 for (i=0; i<size; ++i) {
236 vect[i] = shoot(anEngine); 236 vect[i] = shoot(anEngine);
237 } 237 }
238 } 238 }
239 239
240 void RandGeneral::fireArray( const int size, d 240 void RandGeneral::fireArray( const int size, double* vect )
241 { 241 {
242 int i; 242 int i;
243 243
244 for (i=0; i<size; ++i) { 244 for (i=0; i<size; ++i) {
245 vect[i] = fire(); 245 vect[i] = fire();
246 } 246 }
247 } 247 }
248 248
249 std::ostream & RandGeneral::put ( std::ostream 249 std::ostream & RandGeneral::put ( std::ostream & os ) const {
250 long pr=os.precision(20); << 250 int pr=os.precision(20);
251 std::vector<unsigned long> t(2); 251 std::vector<unsigned long> t(2);
252 os << " " << name() << "\n"; 252 os << " " << name() << "\n";
253 os << "Uvec" << "\n"; 253 os << "Uvec" << "\n";
254 os << nBins << " " << oneOverNbins << " " << 254 os << nBins << " " << oneOverNbins << " " << InterpolationType << "\n";
255 t = DoubConv::dto2longs(oneOverNbins); 255 t = DoubConv::dto2longs(oneOverNbins);
256 os << t[0] << " " << t[1] << "\n"; 256 os << t[0] << " " << t[1] << "\n";
257 assert (static_cast<int>(theIntegralPdf.size 257 assert (static_cast<int>(theIntegralPdf.size())==nBins+1);
258 for (unsigned int i=0; i<theIntegralPdf.size 258 for (unsigned int i=0; i<theIntegralPdf.size(); ++i) {
259 t = DoubConv::dto2longs(theIntegralPdf[i]) 259 t = DoubConv::dto2longs(theIntegralPdf[i]);
260 os << theIntegralPdf[i] << " " << t[0] << 260 os << theIntegralPdf[i] << " " << t[0] << " " << t[1] << "\n";
261 } 261 }
262 os.precision(pr); 262 os.precision(pr);
263 return os; 263 return os;
264 } 264 }
265 265
266 std::istream & RandGeneral::get ( std::istream 266 std::istream & RandGeneral::get ( std::istream & is ) {
267 std::string inName; 267 std::string inName;
268 is >> inName; 268 is >> inName;
269 if (inName != name()) { 269 if (inName != name()) {
270 is.clear(std::ios::badbit | is.rdstate()); 270 is.clear(std::ios::badbit | is.rdstate());
271 std::cerr << "Mismatch when expecting to r 271 std::cerr << "Mismatch when expecting to read state of a "
272 << name() << " distribution\n" 272 << name() << " distribution\n"
273 << "Name found was " << inName 273 << "Name found was " << inName
274 << "\nistream is left in the badbit st 274 << "\nistream is left in the badbit state\n";
275 return is; 275 return is;
276 } 276 }
277 if (possibleKeywordInput(is, "Uvec", nBins)) 277 if (possibleKeywordInput(is, "Uvec", nBins)) {
278 std::vector<unsigned long> t(2); 278 std::vector<unsigned long> t(2);
279 is >> nBins >> oneOverNbins >> Interpolati 279 is >> nBins >> oneOverNbins >> InterpolationType;
280 is >> t[0] >> t[1]; oneOverNbins = DoubCon 280 is >> t[0] >> t[1]; oneOverNbins = DoubConv::longs2double(t);
281 theIntegralPdf.resize(nBins+1); 281 theIntegralPdf.resize(nBins+1);
282 for (unsigned int i=0; i<theIntegralPdf.si 282 for (unsigned int i=0; i<theIntegralPdf.size(); ++i) {
283 is >> theIntegralPdf[i] >> t[0] >> t[1]; 283 is >> theIntegralPdf[i] >> t[0] >> t[1];
284 theIntegralPdf[i] = DoubConv::longs2doub 284 theIntegralPdf[i] = DoubConv::longs2double(t);
285 } 285 }
286 return is; 286 return is;
287 } 287 }
288 // is >> nBins encompassed by possibleKeywor 288 // is >> nBins encompassed by possibleKeywordInput
289 is >> oneOverNbins >> InterpolationType; 289 is >> oneOverNbins >> InterpolationType;
290 theIntegralPdf.resize(nBins+1); 290 theIntegralPdf.resize(nBins+1);
291 for (unsigned int i=0; i<theIntegralPdf.size 291 for (unsigned int i=0; i<theIntegralPdf.size(); ++i) is >> theIntegralPdf[i];
292 return is; 292 return is;
293 } 293 }
294 294
295 } // namespace CLHEP 295 } // namespace CLHEP
296 296