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