Geant4  10.03
RandGeneral.cc
Go to the documentation of this file.
1 // $Id:$
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/RandGeneral.h"
51 #include "CLHEP/Random/DoubConv.h"
52 #include <cassert>
53 
54 namespace CLHEP {
55 
56 std::string RandGeneral::name() const {return "RandGeneral";}
57 HepRandomEngine & RandGeneral::engine() {return *localEngine;}
58 
59 
61 // Constructors
63 
64 RandGeneral::RandGeneral( const double* aProbFunc,
65  int theProbSize,
66  int IntType )
67  : HepRandom(),
68  localEngine(HepRandom::getTheEngine(), do_nothing_deleter()),
69  nBins(theProbSize),
70  InterpolationType(IntType)
71 {
72  prepareTable(aProbFunc);
73 }
74 
75 RandGeneral::RandGeneral(HepRandomEngine& anEngine,
76  const double* aProbFunc,
77  int theProbSize,
78  int IntType )
79 : HepRandom(),
80  localEngine(&anEngine, do_nothing_deleter()),
81  nBins(theProbSize),
82  InterpolationType(IntType)
83 {
84  prepareTable(aProbFunc);
85 }
86 
87 RandGeneral::RandGeneral(HepRandomEngine* anEngine,
88  const double* aProbFunc,
89  int theProbSize,
90  int IntType )
91 : HepRandom(),
92  localEngine(anEngine),
93  nBins(theProbSize),
94  InterpolationType(IntType)
95 {
96  prepareTable(aProbFunc);
97 }
98 
99 void RandGeneral::prepareTable(const double* aProbFunc) {
100 //
101 // Private method called only by constructors. Prepares theIntegralPdf.
102 //
103  if (nBins < 1) {
104  std::cerr <<
105  "RandGeneral constructed with no bins - will use flat distribution\n";
106  useFlatDistribution();
107  return;
108  }
109 
110  theIntegralPdf.resize(nBins+1);
111  theIntegralPdf[0] = 0;
112  int ptn;
113  double weight;
114 
115  for ( ptn = 0; ptn<nBins; ++ptn ) {
116  weight = aProbFunc[ptn];
117  if ( weight < 0 ) {
118  // We can't stomach negative bin contents, they invalidate the
119  // search algorithm when the distribution is fired.
120  std::cerr <<
121  "RandGeneral constructed with negative-weight bin " << ptn <<
122  " = " << weight << " \n -- will substitute 0 weight \n";
123  weight = 0;
124  }
125  // std::cout << ptn << " " << weight << " " << theIntegralPdf[ptn] << "\n";
126  theIntegralPdf[ptn+1] = theIntegralPdf[ptn] + weight;
127  }
128 
129  if ( theIntegralPdf[nBins] <= 0 ) {
130  std::cerr <<
131  "RandGeneral constructed nothing in bins - will use flat distribution\n";
132  useFlatDistribution();
133  return;
134  }
135 
136  for ( ptn = 0; ptn < nBins+1; ++ptn ) {
137  theIntegralPdf[ptn] /= theIntegralPdf[nBins];
138  // std::cout << ptn << " " << theIntegralPdf[ptn] << "\n";
139  }
140 
141  // And another useful variable is ...
142  oneOverNbins = 1.0 / nBins;
143 
144  // One last chore:
145 
146  if ( (InterpolationType != 0) && (InterpolationType != 1) ) {
147  std::cerr <<
148  "RandGeneral does not recognize IntType " << InterpolationType
149  << "\n Will use type 0 (continuous linear interpolation \n";
150  InterpolationType = 0;
151  }
152 
153 } // prepareTable()
154 
155 void RandGeneral::useFlatDistribution() {
156 //
157 // Private method called only by prepareTables in case of user error.
158 //
159  nBins = 1;
160  theIntegralPdf.resize(2);
161  theIntegralPdf[0] = 0;
162  theIntegralPdf[1] = 1;
163  oneOverNbins = 1.0;
164  return;
165 
166 } // UseFlatDistribution()
167 
169 // Destructor
171 
172 RandGeneral::~RandGeneral() {
173 }
174 
175 
177 // mapRandom(rand)
179 
180 double RandGeneral::mapRandom(double rand) const {
181 //
182 // Private method to take the random (however it is created) and map it
183 // according to the distribution.
184 //
185 
186  int nbelow = 0; // largest k such that I[k] is known to be <= rand
187  int nabove = nBins; // largest k such that I[k] is known to be > rand
188  int middle;
189 
190  while (nabove > nbelow+1) {
191  middle = (nabove + nbelow+1)>>1;
192  if (rand >= theIntegralPdf[middle]) {
193  nbelow = middle;
194  } else {
195  nabove = middle;
196  }
197  } // after this loop, nabove is always nbelow+1 and they straddle rad:
198  assert ( nabove == nbelow+1 );
199  assert ( theIntegralPdf[nbelow] <= rand );
200  assert ( theIntegralPdf[nabove] >= rand );
201  // If a defective engine produces rand=1, that will
202  // still give sensible results so we relax the > rand assertion
203 
204  if ( InterpolationType == 1 ) {
205 
206  return nbelow * oneOverNbins;
207 
208  } else {
209 
210  double binMeasure = theIntegralPdf[nabove] - theIntegralPdf[nbelow];
211  // binMeasure is always aProbFunc[nbelow],
212  // but we don't have aProbFunc any more so we subtract.
213 
214  if ( binMeasure == 0 ) {
215  // rand lies right in a bin of measure 0. Simply return the center
216  // of the range of that bin. (Any value between k/N and (k+1)/N is
217  // equally good, in this rare case.)
218  return (nbelow + .5) * oneOverNbins;
219  }
220 
221  double binFraction = (rand - theIntegralPdf[nbelow]) / binMeasure;
222 
223  return (nbelow + binFraction) * oneOverNbins;
224  }
225 
226 } // mapRandom(rand)
227 
228 void RandGeneral::shootArray( HepRandomEngine* anEngine,
229  const int size, double* vect )
230 {
231  int i;
232 
233  for (i=0; i<size; ++i) {
234  vect[i] = shoot(anEngine);
235  }
236 }
237 
238 void RandGeneral::fireArray( const int size, double* vect )
239 {
240  int i;
241 
242  for (i=0; i<size; ++i) {
243  vect[i] = fire();
244  }
245 }
246 
247 std::ostream & RandGeneral::put ( std::ostream & os ) const {
248  int pr=os.precision(20);
249  std::vector<unsigned long> t(2);
250  os << " " << name() << "\n";
251  os << "Uvec" << "\n";
252  os << nBins << " " << oneOverNbins << " " << InterpolationType << "\n";
253  t = DoubConv::dto2longs(oneOverNbins);
254  os << t[0] << " " << t[1] << "\n";
255  assert (static_cast<int>(theIntegralPdf.size())==nBins+1);
256  for (unsigned int i=0; i<theIntegralPdf.size(); ++i) {
257  t = DoubConv::dto2longs(theIntegralPdf[i]);
258  os << theIntegralPdf[i] << " " << t[0] << " " << t[1] << "\n";
259  }
260  os.precision(pr);
261  return os;
262 }
263 
264 std::istream & RandGeneral::get ( std::istream & is ) {
265  std::string inName;
266  is >> inName;
267  if (inName != name()) {
268  is.clear(std::ios::badbit | is.rdstate());
269  std::cerr << "Mismatch when expecting to read state of a "
270  << name() << " distribution\n"
271  << "Name found was " << inName
272  << "\nistream is left in the badbit state\n";
273  return is;
274  }
275  if (possibleKeywordInput(is, "Uvec", nBins)) {
276  std::vector<unsigned long> t(2);
277  is >> nBins >> oneOverNbins >> InterpolationType;
278  is >> t[0] >> t[1]; oneOverNbins = DoubConv::longs2double(t);
279  theIntegralPdf.resize(nBins+1);
280  for (unsigned int i=0; i<theIntegralPdf.size(); ++i) {
281  is >> theIntegralPdf[i] >> t[0] >> t[1];
282  theIntegralPdf[i] = DoubConv::longs2double(t);
283  }
284  return is;
285  }
286  // is >> nBins encompassed by possibleKeywordInput
287  is >> oneOverNbins >> InterpolationType;
288  theIntegralPdf.resize(nBins+1);
289  for (unsigned int i=0; i<theIntegralPdf.size(); ++i) is >> theIntegralPdf[i];
290  return is;
291 }
292 
293 } // namespace CLHEP
ThreeVector shoot(const G4int Ap, const G4int Af)
static int engine(pchar, pchar, double &, pchar &, const dic_type &)
Definition: Evaluator.cc:358
const char * name(G4int ptype)