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 
55 namespace CLHEP {
56 
57 std::string RandGeneral::name() const {return "RandGeneral";}
58 HepRandomEngine & RandGeneral::engine() {return *localEngine;}
59 
60 
62 // Constructors
64 
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 
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 
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  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 
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 
170 // Destructor
172 
174 }
175 
176 
178 // mapRandom(rand)
180 
181 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 
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 
241 void 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 
250 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 
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
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
std::string name() const
Definition: RandGeneral.cc:57
virtual ~RandGeneral()
Definition: RandGeneral.cc:173
static double longs2double(const std::vector< unsigned long > &v)
Definition: DoubConv.cc:106
std::istream & get(std::istream &is)
Definition: RandGeneral.cc:277
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
HepRandomEngine & engine()
Definition: RandGeneral.cc:58
void shootArray(const int size, double *vect)
std::ostream & put(std::ostream &os) const
Definition: RandGeneral.cc:250
void fireArray(const int size, double *vect)
Definition: RandGeneral.cc:241