Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 #if __clang__
30  #if ((defined(G4MULTITHREADED) && !defined(G4USE_STD11)) || \
31  !__has_feature(cxx_thread_local)) || !__has_feature(c_atomic)
32  #define CLANG_NOSTDTLS
33  #endif
34 #endif
35 
36 #if (defined(G4MULTITHREADED) && \
37  (!defined(G4USE_STD11) || (defined(CLANG_NOSTDTLS) || defined(__INTEL_COMPILER))))
38 
39 #include "G4MTRandGeneral.hh"
40 
42 // Constructors
44 
46  G4int theProbSize,
47  G4int IntType )
48  : deleteEngine(false),
49  nBins(theProbSize),
50  InterpolationType(IntType)
51 {
52  localEngine = G4MTHepRandom::getTheEngine();
53  prepareTable(aProbFunc);
54 }
55 
57  const G4double* aProbFunc,
58  G4int theProbSize,
59  G4int IntType )
60 : localEngine(&anEngine),
61  deleteEngine(false),
62  nBins(theProbSize),
63  InterpolationType(IntType)
64 {
65  prepareTable(aProbFunc);
66 }
67 
69  const G4double* aProbFunc,
70  G4int theProbSize,
71  G4int IntType )
72 : localEngine(anEngine),
73  deleteEngine(true),
74  nBins(theProbSize),
75  InterpolationType(IntType)
76 {
77  prepareTable(aProbFunc);
78 }
79 
80 void G4MTRandGeneral::prepareTable(const G4double* aProbFunc)
81 {
82  //
83  // Private method called only by constructors. Prepares theIntegralPdf.
84  //
85  if (nBins < 1) {
86  std::cerr <<
87  "G4MTRandGeneral constructed with no bins - will use flat distribution\n";
88  useFlatDistribution();
89  return;
90  }
91 
92  theIntegralPdf.resize(nBins+1);
93  theIntegralPdf[0] = 0;
94  G4int ptn;
95  G4double weight;
96 
97  for ( ptn = 0; ptn<nBins; ++ptn ) {
98  weight = aProbFunc[ptn];
99  if ( weight < 0 ) {
100  // We can't stomach negative bin contents, they invalidate the
101  // search algorithm when the distribution is fired.
102  std::cerr <<
103  "G4MTRandGeneral constructed with negative-weight bin " << ptn <<
104  " = " << weight << " \n -- will substitute 0 weight \n";
105  weight = 0;
106  }
107  // std::cout << ptn << " " << weight << " " << theIntegralPdf[ptn] << "\n";
108  theIntegralPdf[ptn+1] = theIntegralPdf[ptn] + weight;
109  }
110 
111  if ( theIntegralPdf[nBins] <= 0 ) {
112  std::cerr <<
113  "G4MTRandGeneral constructed nothing in bins - will use flat distribution\n";
114  useFlatDistribution();
115  return;
116  }
117 
118  for ( ptn = 0; ptn < nBins+1; ++ptn ) {
119  theIntegralPdf[ptn] /= theIntegralPdf[nBins];
120  // std::cout << ptn << " " << theIntegralPdf[ptn] << "\n";
121  }
122 
123  // And another useful variable is ...
124  oneOverNbins = 1.0 / nBins;
125 
126  // One last chore:
127 
128  if ( (InterpolationType != 0) && (InterpolationType != 1) ) {
129  std::cerr <<
130  "G4MTRandGeneral does not recognize IntType " << InterpolationType
131  << "\n Will use type 0 (continuous linear interpolation \n";
132  InterpolationType = 0;
133  }
134 
135 } // prepareTable()
136 
137 void G4MTRandGeneral::useFlatDistribution()
138 {
139  //
140  // Private method called only by prepareTables in case of user error.
141  //
142  nBins = 1;
143  theIntegralPdf.resize(2);
144  theIntegralPdf[0] = 0;
145  theIntegralPdf[1] = 1;
146  oneOverNbins = 1.0;
147  return;
148 
149 } // UseFlatDistribution()
150 
152 // Destructor
154 
156 {
157  if ( deleteEngine ) delete localEngine;
158 }
159 
160 
162 // mapRandom(rand)
164 
165 G4double G4MTRandGeneral::mapRandom(G4double rand) const
166 {
167  //
168  // Private method to take the random (however it is created) and map it
169  // according to the distribution.
170  //
171  G4int nbelow = 0; // largest k such that I[k] is known to be <= rand
172  G4int nabove = nBins; // largest k such that I[k] is known to be > rand
173  G4int middle;
174 
175  while (nabove > nbelow+1) {
176  middle = (nabove + nbelow+1)>>1;
177  if (rand >= theIntegralPdf[middle]) {
178  nbelow = middle;
179  } else {
180  nabove = middle;
181  }
182  } // after this loop, nabove is always nbelow+1 and they straddle rad:
183  // assert ( nabove == nbelow+1 );
184  // assert ( theIntegralPdf[nbelow] <= rand );
185  // assert ( theIntegralPdf[nabove] >= rand );
186  // If a defective engine produces rand=1, that will
187  // still give sensible results so we relax the > rand assertion
188 
189  if ( InterpolationType == 1 ) {
190 
191  return nbelow * oneOverNbins;
192 
193  } else {
194 
195  G4double binMeasure = theIntegralPdf[nabove] - theIntegralPdf[nbelow];
196  // binMeasure is always aProbFunc[nbelow],
197  // but we don't have aProbFunc any more so we subtract.
198 
199  if ( binMeasure == 0 ) {
200  // rand lies right in a bin of measure 0. Simply return the center
201  // of the range of that bin. (Any value between k/N and (k+1)/N is
202  // equally good, in this rare case.)
203  return (nbelow + .5) * oneOverNbins;
204  }
205 
206  G4double binFraction = (rand - theIntegralPdf[nbelow]) / binMeasure;
207 
208  return (nbelow + binFraction) * oneOverNbins;
209  }
210 
211 } // mapRandom(rand)
212 
213 
215  const G4int size, G4double* vect )
216 {
217  for (G4int i=0; i<size; ++i)
218  vect[i] = shoot(anEngine);
219 }
220 
221 void G4MTRandGeneral::fireArray( const G4int size, G4double* vect )
222 {
223  for (G4int i=0; i<size; ++i)
224  vect[i] = fire();
225 }
226 #endif
static CLHEP::HepRandomEngine * getTheEngine()
void shootArray(const G4int size, G4double *vect)
G4double fire()
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