| Geant4 Cross Reference |
1 // -*- C++ -*-
2 //
3 // -----------------------------------------------------------------------
4 // HEP Random
5 // --- RandGeneral ---
6 // class implementation file
7 // -----------------------------------------------------------------------
8
9 // =======================================================================
10 // S.Magni & G.Pieri - Created: 5th September 1995
11 // G.Cosmo - Added constructor using default engine from the
12 // static generator. Simplified shoot() and
13 // shootArray() (not needed in principle!): 20th Aug 1998
14 // M.G.Pia & G.Cosmo - Fixed bug in computation of theIntegralPdf in
15 // two constructors: 5th Jan 1999
16 // S.Magni & G.Pieri - Added linear interpolation: 24th Mar 1999
17 // M. Fischler - General cleanup: 14th May 1999
18 // + Eliminated constructor code replication by factoring
19 // common code into prepareTable.
20 // + Eliminated fire/shoot code replication by factoring
21 // out common code into mapRandom.
22 // + A couple of methods are moved inline to avoid a
23 // speed cost for factoring out mapRandom: fire()
24 // and shoot(anEngine).
25 // + Inserted checks for negative weight and zero total
26 // weight in the bins.
27 // + Modified the binary search loop to avoid incorrect
28 // behavior when rand is example on a boundary.
29 // + Moved the check of InterpolationType up into
30 // the constructor. A type other than 0 or 1
31 // will give the interpolated distribution (instead of
32 // a distribution that always returns 0).
33 // + Modified the computation of the returned value
34 // to use algeraic simplification to improve speed.
35 // Eliminated two of the three divisionns, made
36 // use of the fact that nabove-nbelow is always 1, etc.
37 // + Inserted a check for rand hitting the boundary of a
38 // zero-width bin, to avoid dividing 0/0.
39 // M. Fischler - Minor correction in assert 31 July 2001
40 // + changed from assert (above = below+1) to ==
41 // M Fischler - put and get to/from streams 12/15/04
42 // + Modifications to use a vector as theIntegraPdf
43 // M Fischler - put/get to/from streams uses pairs of ulongs when
44 // + storing doubles avoid problems with precision
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 "RandGeneral";}
59 HepRandomEngine & RandGeneral::engine() {return *localEngine;}
60
61
62 //////////////////
63 // Constructors
64 //////////////////
65
66 RandGeneral::RandGeneral( const double* aProbFunc,
67 int theProbSize,
68 int IntType )
69 : HepRandom(),
70 localEngine(HepRandom::getTheEngine(), do_nothing_deleter()),
71 nBins(theProbSize),
72 InterpolationType(IntType)
73 {
74 prepareTable(aProbFunc);
75 }
76
77 RandGeneral::RandGeneral(HepRandomEngine& anEngine,
78 const double* aProbFunc,
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* anEngine,
90 const double* aProbFunc,
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* aProbFunc) {
102 //
103 // Private method called only by constructors. Prepares theIntegralPdf.
104 //
105 if (nBins < 1) {
106 std::cerr <<
107 "RandGeneral constructed with no bins - will use flat distribution\n";
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, they invalidate the
121 // search algorithm when the distribution is fired.
122 std::cerr <<
123 "RandGeneral constructed with negative-weight bin " << ptn <<
124 " = " << weight << " \n -- will substitute 0 weight \n";
125 weight = 0;
126 }
127 // std::cout << ptn << " " << weight << " " << theIntegralPdf[ptn] << "\n";
128 theIntegralPdf[ptn+1] = theIntegralPdf[ptn] + weight;
129 }
130
131 if ( theIntegralPdf[nBins] <= 0 ) {
132 std::cerr <<
133 "RandGeneral constructed nothing in bins - will use flat distribution\n";
134 useFlatDistribution();
135 return;
136 }
137
138 for ( ptn = 0; ptn < nBins+1; ++ptn ) {
139 theIntegralPdf[ptn] /= theIntegralPdf[nBins];
140 // std::cout << ptn << " " << theIntegralPdf[ptn] << "\n";
141 }
142
143 // And another useful variable is ...
144 oneOverNbins = 1.0 / nBins;
145
146 // One last chore:
147
148 if ( (InterpolationType != 0) && (InterpolationType != 1) ) {
149 std::cerr <<
150 "RandGeneral does not recognize IntType " << InterpolationType
151 << "\n Will use type 0 (continuous linear interpolation \n";
152 InterpolationType = 0;
153 }
154
155 } // prepareTable()
156
157 void RandGeneral::useFlatDistribution() {
158 //
159 // Private method called only by prepareTables in case of user error.
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) const {
183 //
184 // Private method to take the random (however it is created) and map it
185 // according to the distribution.
186 //
187
188 int nbelow = 0; // largest k such that I[k] is known to be <= rand
189 int nabove = nBins; // largest k such that I[k] is known to be > rand
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 nbelow+1 and they straddle rad:
200 assert ( nabove == nbelow+1 );
201 assert ( theIntegralPdf[nbelow] <= rand );
202 assert ( theIntegralPdf[nabove] >= rand );
203 // If a defective engine produces rand=1, that will
204 // still give sensible results so we relax the > rand assertion
205
206 if ( InterpolationType == 1 ) {
207
208 return nbelow * oneOverNbins;
209
210 } else {
211
212 double binMeasure = theIntegralPdf[nabove] - theIntegralPdf[nbelow];
213 // binMeasure is always aProbFunc[nbelow],
214 // but we don't have aProbFunc any more so we subtract.
215
216 if ( binMeasure == 0 ) {
217 // rand lies right in a bin of measure 0. Simply return the center
218 // of the range of that bin. (Any value between k/N and (k+1)/N is
219 // equally good, in this rare case.)
220 return (nbelow + .5) * oneOverNbins;
221 }
222
223 double binFraction = (rand - theIntegralPdf[nbelow]) / binMeasure;
224
225 return (nbelow + binFraction) * oneOverNbins;
226 }
227
228 } // mapRandom(rand)
229
230 void RandGeneral::shootArray( HepRandomEngine* anEngine,
231 const int size, double* vect )
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, double* vect )
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 & os ) const {
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 << " " << InterpolationType << "\n";
255 t = DoubConv::dto2longs(oneOverNbins);
256 os << t[0] << " " << t[1] << "\n";
257 assert (static_cast<int>(theIntegralPdf.size())==nBins+1);
258 for (unsigned int i=0; i<theIntegralPdf.size(); ++i) {
259 t = DoubConv::dto2longs(theIntegralPdf[i]);
260 os << theIntegralPdf[i] << " " << t[0] << " " << t[1] << "\n";
261 }
262 os.precision(pr);
263 return os;
264 }
265
266 std::istream & RandGeneral::get ( std::istream & is ) {
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 read state of a "
272 << name() << " distribution\n"
273 << "Name found was " << inName
274 << "\nistream is left in the badbit state\n";
275 return is;
276 }
277 if (possibleKeywordInput(is, "Uvec", nBins)) {
278 std::vector<unsigned long> t(2);
279 is >> nBins >> oneOverNbins >> InterpolationType;
280 is >> t[0] >> t[1]; oneOverNbins = DoubConv::longs2double(t);
281 theIntegralPdf.resize(nBins+1);
282 for (unsigned int i=0; i<theIntegralPdf.size(); ++i) {
283 is >> theIntegralPdf[i] >> t[0] >> t[1];
284 theIntegralPdf[i] = DoubConv::longs2double(t);
285 }
286 return is;
287 }
288 // is >> nBins encompassed by possibleKeywordInput
289 is >> oneOverNbins >> InterpolationType;
290 theIntegralPdf.resize(nBins+1);
291 for (unsigned int i=0; i<theIntegralPdf.size(); ++i) is >> theIntegralPdf[i];
292 return is;
293 }
294
295 } // namespace CLHEP
296