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