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