Geant4_10
G4MTRandGeneral.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id:$
28 //
29 #ifdef G4MULTITHREADED
30 
31 #include "G4MTRandGeneral.hh"
32 
34 // Constructors
36 
38  G4int theProbSize,
39  G4int IntType )
40  : deleteEngine(false),
41  nBins(theProbSize),
42  InterpolationType(IntType)
43 {
44  localEngine = G4MTHepRandom::getTheEngine();
45  prepareTable(aProbFunc);
46 }
47 
49  const G4double* aProbFunc,
50  G4int theProbSize,
51  G4int IntType )
52 : localEngine(&anEngine),
53  deleteEngine(false),
54  nBins(theProbSize),
55  InterpolationType(IntType)
56 {
57  prepareTable(aProbFunc);
58 }
59 
61  const G4double* aProbFunc,
62  G4int theProbSize,
63  G4int IntType )
64 : localEngine(anEngine),
65  deleteEngine(true),
66  nBins(theProbSize),
67  InterpolationType(IntType)
68 {
69  prepareTable(aProbFunc);
70 }
71 
72 void G4MTRandGeneral::prepareTable(const G4double* aProbFunc)
73 {
74  //
75  // Private method called only by constructors. Prepares theIntegralPdf.
76  //
77  if (nBins < 1) {
78  std::cerr <<
79  "G4MTRandGeneral constructed with no bins - will use flat distribution\n";
80  useFlatDistribution();
81  return;
82  }
83 
84  theIntegralPdf.resize(nBins+1);
85  theIntegralPdf[0] = 0;
86  register G4int ptn;
87  register G4double weight;
88 
89  for ( ptn = 0; ptn<nBins; ++ptn ) {
90  weight = aProbFunc[ptn];
91  if ( weight < 0 ) {
92  // We can't stomach negative bin contents, they invalidate the
93  // search algorithm when the distribution is fired.
94  std::cerr <<
95  "G4MTRandGeneral constructed with negative-weight bin " << ptn <<
96  " = " << weight << " \n -- will substitute 0 weight \n";
97  weight = 0;
98  }
99  // std::cout << ptn << " " << weight << " " << theIntegralPdf[ptn] << "\n";
100  theIntegralPdf[ptn+1] = theIntegralPdf[ptn] + weight;
101  }
102 
103  if ( theIntegralPdf[nBins] <= 0 ) {
104  std::cerr <<
105  "G4MTRandGeneral constructed nothing in bins - will use flat distribution\n";
106  useFlatDistribution();
107  return;
108  }
109 
110  for ( ptn = 0; ptn < nBins+1; ++ptn ) {
111  theIntegralPdf[ptn] /= theIntegralPdf[nBins];
112  // std::cout << ptn << " " << theIntegralPdf[ptn] << "\n";
113  }
114 
115  // And another useful variable is ...
116  oneOverNbins = 1.0 / nBins;
117 
118  // One last chore:
119 
120  if ( (InterpolationType != 0) && (InterpolationType != 1) ) {
121  std::cerr <<
122  "G4MTRandGeneral does not recognize IntType " << InterpolationType
123  << "\n Will use type 0 (continuous linear interpolation \n";
124  InterpolationType = 0;
125  }
126 
127 } // prepareTable()
128 
129 void G4MTRandGeneral::useFlatDistribution()
130 {
131  //
132  // Private method called only by prepareTables in case of user error.
133  //
134  nBins = 1;
135  theIntegralPdf.resize(2);
136  theIntegralPdf[0] = 0;
137  theIntegralPdf[1] = 1;
138  oneOverNbins = 1.0;
139  return;
140 
141 } // UseFlatDistribution()
142 
144 // Destructor
146 
148 {
149  if ( deleteEngine ) delete localEngine;
150 }
151 
152 
154 // mapRandom(rand)
156 
157 G4double G4MTRandGeneral::mapRandom(G4double rand) const
158 {
159  //
160  // Private method to take the random (however it is created) and map it
161  // according to the distribution.
162  //
163  G4int nbelow = 0; // largest k such that I[k] is known to be <= rand
164  G4int nabove = nBins; // largest k such that I[k] is known to be > rand
165  G4int middle;
166 
167  while (nabove > nbelow+1) {
168  middle = (nabove + nbelow+1)>>1;
169  if (rand >= theIntegralPdf[middle]) {
170  nbelow = middle;
171  } else {
172  nabove = middle;
173  }
174  } // after this loop, nabove is always nbelow+1 and they straddle rad:
175  // assert ( nabove == nbelow+1 );
176  // assert ( theIntegralPdf[nbelow] <= rand );
177  // assert ( theIntegralPdf[nabove] >= rand );
178  // If a defective engine produces rand=1, that will
179  // still give sensible results so we relax the > rand assertion
180 
181  if ( InterpolationType == 1 ) {
182 
183  return nbelow * oneOverNbins;
184 
185  } else {
186 
187  G4double binMeasure = theIntegralPdf[nabove] - theIntegralPdf[nbelow];
188  // binMeasure is always aProbFunc[nbelow],
189  // but we don't have aProbFunc any more so we subtract.
190 
191  if ( binMeasure == 0 ) {
192  // rand lies right in a bin of measure 0. Simply return the center
193  // of the range of that bin. (Any value between k/N and (k+1)/N is
194  // equally good, in this rare case.)
195  return (nbelow + .5) * oneOverNbins;
196  }
197 
198  G4double binFraction = (rand - theIntegralPdf[nbelow]) / binMeasure;
199 
200  return (nbelow + binFraction) * oneOverNbins;
201  }
202 
203 } // mapRandom(rand)
204 
205 
207  const G4int size, G4double* vect )
208 {
209  for (register G4int i=0; i<size; ++i)
210  vect[i] = shoot(anEngine);
211 }
212 
213 void G4MTRandGeneral::fireArray( const G4int size, G4double* vect )
214 {
215  for (register G4int i=0; i<size; ++i)
216  vect[i] = fire();
217 }
218 #endif
static CLHEP::HepRandomEngine * getTheEngine()
void shootArray(const G4int size, G4double *vect)
G4double fire()
double weight
Definition: plottest35.C:25
int G4int
Definition: G4Types.hh:78
void fireArray(const G4int size, G4double *vect)
G4MTRandGeneral(const G4double *aProbFunc, G4int theProbSize, G4int IntType=0)
virtual ~G4MTRandGeneral()
G4double shoot()
double G4double
Definition: G4Types.hh:76