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
CLHEP::HepRandomEngine
Definition: Matrix/CLHEP/Random/RandomEngine.h:55
CLHEP::RandGeneral::fireArray
void fireArray(const int size, double *vect)
Definition: RandGeneral.cc:241
CLHEP::RandGeneral::~RandGeneral
virtual ~RandGeneral()
Definition: RandGeneral.cc:173
is
HepRotation and so forth isNear() norm2() rectify() static Rotation row1 row4(To avoid bloat in the code pulled in for programs which don 't use all these features, we split the implementation .cc files. Only isNear() goes into the original Rotation.cc) --------------------------------------- HepAxisAngle and HepEulerAngles classes --------------------------------------- These classes are very useful and simple structures for holding the result of a nice intuituve decomposition of a rotation there is no longer much content in the distinct ZOOM PhysicsVectors library The only content left in the library is the object files representing the various Exception objects When we build the CLHEP classes for the ZOOM we will set up so as to use ZOOM SpaceVector is(but we can disable namespace usage and most of our users do so at this point). What I do is leave Hep3Vector in the global namespace
CLHEP::DoubConv::longs2double
static double longs2double(const std::vector< unsigned long > &v)
Definition: DoubConv.cc:106
CLHEP::RandGeneral::RandGeneral
RandGeneral(const double *aProbFunc, int theProbSize, int IntType=0)
Definition: RandGeneral.cc:65
size
user code seldom needs to call this function directly ZMerrno whether or not they are still recorded ZMerrno size() Return the(integer) number of ZMthrow 'n exceptions currently recorded. 5) ZMerrno.clear() Set an internal counter to zero. This counter is available(see next function) to user code to track ZMthrow 'n exceptions that have occurred during any arbitrary time interval. 6) ZMerrno.countSinceCleared() Return the(integer) number of ZMthrow 'n exceptions that have been recorded via ZMerrno.write()
CLHEP::RandGeneral::shoot
double shoot()
CLHEP
Definition: ClhepVersion.h:13
CLHEP::do_nothing_deleter
Definition: Matrix/CLHEP/Utility/memory.h:1477
CLHEP::RandGeneral::engine
HepRandomEngine & engine()
Definition: RandGeneral.cc:58
CLHEP::RandGeneral::get
std::istream & get(std::istream &is)
Definition: RandGeneral.cc:277
CLHEP::possibleKeywordInput
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
Definition: Matrix/CLHEP/Random/RandomEngine.h:168
CLHEP::RandGeneral::name
std::string name() const
Definition: RandGeneral.cc:57
i
long i
Definition: JamesRandomSeeding.txt:27
CLHEP::DoubConv::dto2longs
static std::vector< unsigned long > dto2longs(double d)
Definition: DoubConv.cc:90
CLHEP::RandGeneral::fire
double fire()
CLHEP::RandGeneral::shootArray
void shootArray(const int size, double *vect)
CLHEP::HepRandom
Definition: Matrix/CLHEP/Random/Random.h:50
CLHEP::RandGeneral::put
std::ostream & put(std::ostream &os) const
Definition: RandGeneral.cc:250