CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

RandGeneral.cc
Go to the documentation of this file.
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
55namespace CLHEP {
56
57std::string RandGeneral::name() const {return "RandGeneral";}
58HepRandomEngine & RandGeneral::engine() {return *localEngine;}
59
60
62// Constructors
64
65RandGeneral::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
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
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
100void 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 register int ptn;
114 register 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
156void 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
170// Destructor
172
174}
175
176
178// mapRandom(rand)
180
181double 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
232 const int size, double* vect )
233{
234 register int i;
235
236 for (i=0; i<size; ++i) {
237 vect[i] = shoot(anEngine);
238 }
239}
240
241void RandGeneral::fireArray( const int size, double* vect )
242{
243 register int i;
244
245 for (i=0; i<size; ++i) {
246 vect[i] = fire();
247 }
248}
249
250std::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
277std::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
static double longs2double(const std::vector< unsigned long > &v)
Definition: DoubConv.cc:106
static std::vector< unsigned long > dto2longs(double d)
Definition: DoubConv.cc:90
RandGeneral(const double *aProbFunc, int theProbSize, int IntType=0)
Definition: RandGeneral.cc:65
std::istream & get(std::istream &is)
Definition: RandGeneral.cc:277
virtual ~RandGeneral()
Definition: RandGeneral.cc:173
std::ostream & put(std::ostream &os) const
Definition: RandGeneral.cc:250
void fireArray(const int size, double *vect)
Definition: RandGeneral.cc:241
std::string name() const
Definition: RandGeneral.cc:57
HepRandomEngine & engine()
Definition: RandGeneral.cc:58
void shootArray(const int size, double *vect)
bool possibleKeywordInput(IS &is, const std::string &key, T &t)