| Geant4 Cross Reference |
1 // -*- C++ -*- 1
2 //
3 // -------------------------------------------
4 // HEP Random
5 // --- RandLandau ---
6 // class implementation f
7 // -------------------------------------------
8
9 // ===========================================
10 // M Fischler - Created 1/6/2000.
11 //
12 // The key transform() method uses the
13 // This is because I trust that RANLAN
14 // I trust the Bukin-Grozina inverseLan
15 // claimed to be better than 1% accurat
16 //
17 // M Fischler - put and get to/from strea
18 // ===========================================
19
20 #include "CLHEP/Random/RandLandau.h"
21 #include <iostream>
22 #include <cmath> // for std::log()
23
24 namespace CLHEP {
25
26 std::string RandLandau::name() const {return "
27 HepRandomEngine & RandLandau::engine() {return
28
29 RandLandau::~RandLandau() {
30 }
31
32 void RandLandau::shootArray( const int size, d
33
34 {
35 for( double* v = vect; v != vect + size; ++v
36 *v = shoot();
37 }
38
39 void RandLandau::shootArray( HepRandomEngine*
40 const int size, do
41 {
42 for( double* v = vect; v != vect + size; ++v
43 *v = shoot(anEngine);
44 }
45
46 void RandLandau::fireArray( const int size, do
47 {
48 for( double* v = vect; v != vect + size; ++v
49 *v = fire();
50 }
51
52 //
53 // Table of values of inverse Landau, from r =
54 //
55
56 // Since all these are this is static to this
57 // info is establised a priori and not at each
58
59 static const float TABLE_INTERVAL = .001f;
60 static const int TABLE_END = 982;
61 static const float TABLE_MULTIPLIER = 1.0f/TAB
62
63 // Here comes the big (4K bytes) table ---
64 //
65 // inverseLandau[ n ] = the inverse cdf at r =
66 //
67 // Credit CERNLIB for these computations.
68 //
69 // This data is float because the algortihm do
70 // data accuracy. The numbers below .006 or a
71 // non-table-based methods are used below r=.0
72
73 static const float inverseLandau [TABLE_END+1]
74
75 0.0f, // .000
76 0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
77 -2.244733f, -2.204365f, -2.168163f, -2.135219f
78 -2.076740f, -2.050397f, -2.025605f, -2.002150f
79 -1.958612f, -1.938275f, -1.918760f, -1.899984f
80 -1.864385f, -1.847451f, -1.831030f, -1.815083f
81 -1.784473f, -1.769751f, -1.755383f, -1.741346f
82 -1.714187f, -1.701029f, -1.688130f, -1.675477f
83 -1.650858f, -1.638868f, -1.627078f, -1.615477f
84 -1.592811f, -1.581729f, -1.570806f, -1.560034f
85 -1.538919f, -1.528565f, -1.518339f, -1.508237f
86 -1.488386f, -1.478628f, -1.468976f, -1.459428f
87 -1.440626f, -1.431365f, -1.422195f, -1.413111f
88 -1.395194f, -1.386356f, -1.377594f, -1.368906f
89 -1.351746f, -1.343269f, -1.334859f, -1.326512f
90 -1.310006f, -1.301843f, -1.293737f, -1.285688f
91 -1.269752f, -1.261863f, -1.254024f, -1.246235f
92 -1.230800f, -1.223153f, -1.215550f, -1.207990f
93 -1.192999f, -1.185566f, -1.178172f, -1.170817f
94 -1.156220f, -1.148977f, -1.141770f, -1.134598f
95 -1.120354f, -1.113282f, -1.106242f, -1.099233f
96
97 -1.085306f, -1.078388f, -1.071498f, -1.064636f
98 -1.050996f, -1.044215f, -1.037461f, -1.030733f
99 -1.017350f, -1.010695f, -1.004064f, -.997456f,
100 -.984308f, -.977767f, -.971247f, -.964749f, -.
101 -.951813f, -.945375f, -.938957f, -.932558f, -.
102 -.919816f, -.913472f, -.907146f, -.900838f, -.
103 -.888272f, -.882014f, -.875773f, -.869547f, -.
104 -.857142f, -.850963f, -.844798f, -.838648f, -.
105 -.826390f, -.820282f, -.814187f, -.808106f, -.
106 -.795982f, -.789940f, -.783909f, -.777891f, -.
107 -.765889f, -.759906f, -.753934f, -.747973f, -.
108 -.736084f, -.730155f, -.724237f, -.718328f, -.
109 -.706541f, -.700661f, -.694791f, -.688931f, -.
110 -.677236f, -.671402f, -.665576f, -.659759f, -.
111 -.648149f, -.642356f, -.636570f, -.630793f, -.
112 -.619259f, -.613503f, -.607754f, -.602012f, -.
113 -.590548f, -.584825f, -.579109f, -.573399f, -.
114 -.561997f, -.556305f, -.550618f, -.544937f, -.
115 -.533592f, -.527926f, -.522266f, -.516611f, -.
116 -.505315f, -.499674f, -.494037f, -.488405f, -.
117
118 -.477153f, -.471533f, -.465917f, -.460305f, -.
119 -.449092f, -.443491f, -.437893f, -.432299f, -.
120 -.421119f, -.415534f, -.409951f, -.404372f, -.
121 -.393221f, -.387649f, -.382080f, -.376513f, -.
122 -.365387f, -.359826f, -.354268f, -.348712f, -.
123 -.337604f, -.332053f, -.326503f, -.320955f, -.
124 -.309863f, -.304318f, -.298775f, -.293233f, -.
125 -.282152f, -.276613f, -.271074f, -.265536f, -.
126 -.254462f, -.248926f, -.243389f, -.237854f, -.
127 -.226783f, -.221247f, -.215712f, -.210176f, -.
128 -.199105f, -.193568f, -.188032f, -.182495f, -.
129 -.171419f, -.165880f, -.160341f, -.154800f, -.
130 -.143717f, -.138173f, -.132629f, -.127083f, -.
131 -.115989f, -.110439f, -.104889f, -.099336f, -.
132 -.088227f, -.082670f, -.077111f, -.071550f, -.
133 -.060423f, -.054856f, -.049288f, -.043717f, -.
134 -.032569f, -.026991f, -.021411f, -.015828f, -.
135 -.004656f, .000934f, .006527f, .012123f, .
136 .023323f, .028928f, .034535f, .040146f, .04
137 .051376f, .056997f, .062620f, .068247f, .07
138
139 .079511f, .085149f, .090790f, .096435f, .1
140 .107736f, .113392f, .119052f, .124716f, .1
141 .136057f, .141734f, .147414f, .153100f, .1
142 .164483f, .170181f, .175884f, .181592f, .1
143 .193021f, .198743f, .204469f, .210201f, .2
144 .221678f, .227425f, .233177f, .238933f, .2
145 .250463f, .256236f, .262014f, .267798f, .2
146 .279382f, .285183f, .290989f, .296801f, .3
147 .308443f, .314273f, .320109f, .325951f, .3
148 .337654f, .343515f, .349382f, .355255f, .3
149 .367022f, .372915f, .378815f, .384721f, .3
150 .396554f, .402481f, .408415f, .414356f, .4
151 .426260f, .432222f, .438192f, .444169f, .4
152 .456145f, .462144f, .468151f, .474166f, .4
153 .486218f, .492256f, .498302f, .504356f, .5
154 .516488f, .522566f, .528653f, .534747f, .5
155 .546962f, .553082f, .559210f, .565347f, .5
156 .577648f, .583811f, .589983f, .596164f, .6
157 .608554f, .614762f, .620980f, .627207f, .6
158 .639689f, .645945f, .652210f, .658484f, .6
159
160 .671062f, .677366f, .683680f, .690004f, .6
161 .702682f, .709036f, .715400f, .721775f, .7
162 .734556f, .740963f, .747379f, .753807f, .7
163 .766695f, .773155f, .779627f, .786109f, .7
164 .799107f, .805624f, .812151f, .818690f, .8
165 .831803f, .838377f, .844962f, .851560f, .8
166 .864791f, .871425f, .878071f, .884729f, .8
167 .898082f, .904778f, .911486f, .918206f, .9
168 .931686f, .938446f, .945218f, .952003f, .9
169 .965614f, .972439f, .979278f, .986130f, .9
170 .999875f, 1.006769f, 1.013676f, 1.020597f, 1.
171 1.034482f, 1.041446f, 1.048424f, 1.055417f, 1.
172 1.069446f, 1.076482f, 1.083534f, 1.090600f, 1.
173 1.104778f, 1.111889f, 1.119016f, 1.126159f, 1.
174 1.140490f, 1.147679f, 1.154884f, 1.162105f, 1.
175 1.176595f, 1.183864f, 1.191149f, 1.198451f, 1.
176 1.213105f, 1.220457f, 1.227826f, 1.235211f, 1.
177 1.250034f, 1.257471f, 1.264926f, 1.272398f, 1.
178 1.287395f, 1.294921f, 1.302464f, 1.310026f, 1.
179 1.325203f, 1.332819f, 1.340454f, 1.348108f, 1.
180
181 1.363472f, 1.371182f, 1.378912f, 1.386660f, 1.
182 1.402216f, 1.410024f, 1.417851f, 1.425698f, 1.
183 1.441453f, 1.449360f, 1.457288f, 1.465237f, 1.
184 1.481196f, 1.489208f, 1.497240f, 1.505293f, 1.
185 1.521465f, 1.529583f, 1.537723f, 1.545885f, 1.
186 1.562275f, 1.570503f, 1.578754f, 1.587028f, 1.
187 1.603644f, 1.611987f, 1.620353f, 1.628743f, 1.
188 1.645593f, 1.654053f, 1.662538f, 1.671047f, 1.
189 1.688139f, 1.696721f, 1.705329f, 1.713961f, 1.
190 1.731303f, 1.740011f, 1.748746f, 1.757506f, 1.
191 1.775106f, 1.783945f, 1.792810f, 1.801703f, 1.
192 1.819569f, 1.828543f, 1.837545f, 1.846574f, 1.
193 1.864717f, 1.873830f, 1.882972f, 1.892143f, 1.
194 1.910572f, 1.919830f, 1.929117f, 1.938434f, 1.
195 1.957158f, 1.966566f, 1.976004f, 1.985473f, 1.
196 2.004503f, 2.014065f, 2.023659f, 2.033285f, 2.
197 2.052633f, 2.062355f, 2.072110f, 2.081899f, 2.
198 2.101575f, 2.111464f, 2.121386f, 2.131343f, 2.
199 2.151360f, 2.161421f, 2.171517f, 2.181648f, 2.
200 2.202018f, 2.212257f, 2.222533f, 2.232845f, 2.
201
202 2.253582f, 2.264006f, 2.274468f, 2.284968f, 2.
203 2.306084f, 2.316701f, 2.327356f, 2.338051f, 2.
204 2.359562f, 2.370377f, 2.381234f, 2.392131f, 2.
205 2.414051f, 2.425073f, 2.436138f, 2.447246f, 2.
206 2.469591f, 2.480828f, 2.492110f, 2.503436f, 2.
207 2.526222f, 2.537684f, 2.549190f, 2.560743f, 2.
208 2.583989f, 2.595682f, 2.607423f, 2.619212f, 2.
209 2.642936f, 2.654871f, 2.666855f, 2.678890f, 2.
210 2.703110f, 2.715297f, 2.727535f, 2.739825f, 2.
211 2.764563f, 2.777012f, 2.789514f, 2.802070f, 2.
212 2.827347f, 2.840069f, 2.852846f, 2.865680f, 2.
213 2.891518f, 2.904524f, 2.917588f, 2.930712f, 2.
214 2.957136f, 2.970439f, 2.983802f, 2.997227f, 3.
215 3.024263f, 3.037875f, 3.051551f, 3.065290f, 3.
216 3.092965f, 3.106900f, 3.120902f, 3.134971f, 3.
217 3.163312f, 3.177585f, 3.191928f, 3.206340f, 3.
218 3.235378f, 3.250005f, 3.264704f, 3.279477f, 3.
219 3.309244f, 3.324240f, 3.339312f, 3.354461f, 3.
220 3.384992f, 3.400375f, 3.415838f, 3.431381f, 3.
221 3.462711f, 3.478500f, 3.494372f, 3.510328f, 3.
222
223 3.542497f, 3.558711f, 3.575012f, 3.591402f, 3.
224 3.624450f, 3.641111f, 3.657863f, 3.674708f, 3.
225 3.708680f, 3.725809f, 3.743034f, 3.760357f, 3.
226 3.795300f, 3.812921f, 3.830645f, 3.848470f, 3.
227 3.884434f, 3.902574f, 3.920821f, 3.939176f, 3.
228 3.976215f, 3.994901f, 4.013699f, 4.032612f, 4.
229 4.070783f, 4.090045f, 4.109425f, 4.128925f, 4.
230 4.168292f, 4.188160f, 4.208154f, 4.228275f, 4.
231 4.268903f, 4.289413f, 4.310056f, 4.330832f, 4.
232 4.372794f, 4.393982f, 4.415310f, 4.436781f, 4.
233 4.480154f, 4.502060f, 4.524114f, 4.546319f, 4.
234 4.591187f, 4.613854f, 4.636678f, 4.659662f, 4.
235 4.706116f, 4.729590f, 4.753231f, 4.777041f, 4.
236 4.825179f, 4.849511f, 4.874020f, 4.898710f, 4.
237 4.948639f, 4.973883f, 4.999316f, 5.024942f, 5.
238 5.076778f, 5.102993f, 5.129411f, 5.156034f, 5.
239 5.209903f, 5.237156f, 5.264625f, 5.292312f, 5.
240 5.348354f, 5.376714f, 5.405306f, 5.434131f, 5.
241 5.492496f, 5.522042f, 5.551836f, 5.581880f, 5.
242 5.642734f, 5.673552f, 5.704634f, 5.735986f, 5.
243
244 5.799512f, 5.831694f, 5.864161f, 5.896918f, 5.
245 5.963316f, 5.996967f, 6.030925f, 6.065194f, 6.
246 6.134687f, 6.169921f, 6.205486f, 6.241387f, 6.
247 6.314220f, 6.351163f, 6.388465f, 6.426130f, 6.
248 6.502578f, 6.541371f, 6.580553f, 6.620130f, 6.
249 6.700495f, 6.741297f, 6.782520f, 6.824173f, 6.
250 6.908795f, 6.951780f, 6.995225f, 7.039137f, 7.
251 7.128398f, 7.173764f, 7.219632f, 7.266011f, 7.
252 7.360339f, 7.408308f, 7.456827f, 7.505905f, 7.
253 7.605785f, 7.656608f, 7.708035f, 7.760077f, 7.
254 7.866057f, 7.920019f, 7.974647f, 8.029953f, 8.
255 8.142657f, 8.200083f, 8.258245f, 8.317158f, 8.
256 8.437300f, 8.498562f, 8.560641f, 8.623554f, 8.
257 8.751955f, 8.817481f, 8.883916f, 8.951282f, 9.
258 9.088889f, 9.159174f, 9.230477f, 9.302822f, 9.
259 9.450735f, 9.526355f, 9.603118f, 9.681054f, 9.
260 9.840558f, 9.922186f, 10.005107f, 10.089353f
261 10.261958f, 10.350389f, 10.440287f, 10.531693f
262 10.719188f, 10.815362f, 10.913214f, 11.012789f
263 11.217307f, 11.322352f, 11.429325f, 11.538283f
264
265 11.762390f, 11.877664f, 11.995170f, 12.114979f
266 12.361791f, 12.488946f, 12.618708f, 12.751161f
267 13.024498f, 13.165570f, 13.309711f, 13.457026f
268 13.761625f, 13.919145f, 14.080314f, 14.245263f
269 14.587072f, 14.764233f, 14.945778f, 15.131877f
270 15.518470f, 15.719353f, 15.925570f, 16.137345f
271 16.578520f, 16.808433f, 17.044929f, 17.288305f
272 17.796967f, 18.062943f, 18.337176f, 18.620068f
273 19.213574f, 19.525133f, 19.847249f, 20.180480f
274 20.882738f, 21.253102f, 21.637266f, 22.036036f
275 22.880933f, 23.329017f, 23.795634f, 24.281981f
276 25.319207f, 25.873062f, 26.452634f, 27.059789f
277 28.365274f, 29.068370f, 29.808638f, 30.589157f
278 32.285060f, 33.208568f, 34.188705f, 35.230920f
279 37.527131f, 38.796172f, 40.157721f, 41.622399f
280 44.912465f, 46.769077f, 48.792279f, 51.005773f
281 56.123356f, 59.103894f, // .982
282
283 }; // End of the inverseLandau table
284
285 double RandLandau::transform (double r) {
286
287 double u = r * TABLE_MULTIPLIER;
288 int index = int(u);
289 double du = u - index;
290
291 // du is scaled such that the we dont have t
292 // when interpolating.
293
294 // Five cases:
295 // A) Between .070 and .800 the function is
296 // linear interpolation is adequate.
297 // B) Between .007 and .070, and between .80
298 // interpolation is used. This requires
299 // a cubic spline (thus we need .006 and .9
300 // the quadratic interpolation is accurate
301 // C) Below .007 an asymptotic expansion for
302 // (involving two logs) is used; there is
303 // factor.
304 // D) Above .980, a simple pade approximatio
305 // 1/(1-r)), but...
306 // E) the coefficients in that pade are diff
307
308 if ( index >= 70 && index <= 800 ) { // (
309
310 double f0 = inverseLandau [index];
311 double f1 = inverseLandau [index+1];
312 return f0 + du * (f1 - f0);
313
314 } else if ( index >= 7 && index <= 980 ) {
315
316 double f_1 = inverseLandau [index-1];
317 double f0 = inverseLandau [index];
318 double f1 = inverseLandau [index+1];
319 double f2 = inverseLandau [index+2];
320
321 return f0 + du * (f1 - f0 - .25*(1-du)* (f
322
323 } else if ( index < 7 ) { // (C)
324
325 const double n0 = 0.99858950;
326 const double n1 = 34.5213058; const double
327 const double n2 = 17.0854528; const double
328
329 double logr = std::log(r);
330 double x = 1/logr;
331 double x2 = x*x;
332
333 double pade = (n0 + n1*x + n2*x2) / (1.0 +
334
335 return ( - std::log ( -.91893853 - logr )
336
337 } else if ( index <= 999 ) { // (D)
338
339 const double n0 = 1.00060006;
340 const double n1 = 263.991156; const doub
341 const double n2 = 4373.20068; const double
342
343 double x = 1-r;
344 double x2 = x*x;
345
346 return (n0 + n1*x + n2*x2) / (x * (1.0 + d
347
348 } else { // (E)
349
350 const double n0 = 1.00001538;
351 const double n1 = 6075.14119; const doub
352 const double n2 = 734266.409; const double
353
354 double x = 1-r;
355 double x2 = x*x;
356
357 return (n0 + n1*x + n2*x2) / (x * (1.0 + d
358
359 }
360
361 } // transform()
362
363 std::ostream & RandLandau::put ( std::ostream
364 long pr=os.precision(20);
365 os << " " << name() << "\n";
366 os.precision(pr);
367 return os;
368 }
369
370 std::istream & RandLandau::get ( std::istream
371 std::string inName;
372 is >> inName;
373 if (inName != name()) {
374 is.clear(std::ios::badbit | is.rdstate());
375 std::cerr << "Mismatch when expecting to r
376 << name() << " distribution\n"
377 << "Name found was " << inName
378 << "\nistream is left in the badbit st
379 return is;
380 }
381 return is;
382 }
383
384 } // namespace CLHEP
385