CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

RandLandau.cc
Go to the documentation of this file.
1// $Id: RandLandau.cc,v 1.5 2010/06/16 17:24:53 garren Exp $
2// -*- C++ -*-
3//
4// -----------------------------------------------------------------------
5// HEP Random
6// --- RandLandau ---
7// class implementation file
8// -----------------------------------------------------------------------
9
10// =======================================================================
11// M Fischler - Created 1/6/2000.
12//
13// The key transform() method uses the algorithm in CERNLIB.
14// This is because I trust that RANLAN routine more than
15// I trust the Bukin-Grozina inverseLandau, which is not
16// claimed to be better than 1% accurate.
17//
18// M Fischler - put and get to/from streams 12/13/04
19// =======================================================================
20
21#include "CLHEP/Random/defs.h"
23#include <iostream>
24#include <cmath> // for std::log()
25
26namespace CLHEP {
27
28std::string RandLandau::name() const {return "RandLandau";}
29HepRandomEngine & RandLandau::engine() {return *localEngine;}
30
33
34void RandLandau::shootArray( const int size, double* vect )
35
36{
37 for( double* v = vect; v != vect + size; ++v )
38 *v = shoot();
39}
40
42 const int size, double* vect )
43{
44 for( double* v = vect; v != vect + size; ++v )
45 *v = shoot(anEngine);
46}
47
48void RandLandau::fireArray( const int size, double* vect)
49{
50 for( double* v = vect; v != vect + size; ++v )
51 *v = fire();
52}
53
54
55//
56// Table of values of inverse Landau, from r = .060 to .982
57//
58
59// Since all these are this is static to this compilation unit only, the
60// info is establised a priori and not at each invocation.
61
62static const float TABLE_INTERVAL = .001f;
63static const int TABLE_END = 982;
64static const float TABLE_MULTIPLIER = 1.0f/TABLE_INTERVAL;
65
66// Here comes the big (4K bytes) table ---
67//
68// inverseLandau[ n ] = the inverse cdf at r = n*TABLE_INTERVAL = n/1000.
69//
70// Credit CERNLIB for these computations.
71//
72// This data is float because the algortihm does not benefit from further
73// data accuracy. The numbers below .006 or above .982 are moot since
74// non-table-based methods are used below r=.007 and above .980.
75
76static const float inverseLandau [TABLE_END+1] = {
77
780.0f, // .000
790.0f, 0.0f, 0.0f, 0.0f, 0.0f, // .001 - .005
80-2.244733f, -2.204365f, -2.168163f, -2.135219f, -2.104898f, // .006 - .010
81-2.076740f, -2.050397f, -2.025605f, -2.002150f, -1.979866f,
82-1.958612f, -1.938275f, -1.918760f, -1.899984f, -1.881879f, // .020
83-1.864385f, -1.847451f, -1.831030f, -1.815083f, -1.799574f,
84-1.784473f, -1.769751f, -1.755383f, -1.741346f, -1.727620f, // .030
85-1.714187f, -1.701029f, -1.688130f, -1.675477f, -1.663057f,
86-1.650858f, -1.638868f, -1.627078f, -1.615477f, -1.604058f, // .040
87-1.592811f, -1.581729f, -1.570806f, -1.560034f, -1.549407f,
88-1.538919f, -1.528565f, -1.518339f, -1.508237f, -1.498254f, // .050
89-1.488386f, -1.478628f, -1.468976f, -1.459428f, -1.449979f,
90-1.440626f, -1.431365f, -1.422195f, -1.413111f, -1.404112f, // .060
91-1.395194f, -1.386356f, -1.377594f, -1.368906f, -1.360291f,
92-1.351746f, -1.343269f, -1.334859f, -1.326512f, -1.318229f, // .070
93-1.310006f, -1.301843f, -1.293737f, -1.285688f, -1.277693f,
94-1.269752f, -1.261863f, -1.254024f, -1.246235f, -1.238494f, // .080
95-1.230800f, -1.223153f, -1.215550f, -1.207990f, -1.200474f,
96-1.192999f, -1.185566f, -1.178172f, -1.170817f, -1.163500f, // .090
97-1.156220f, -1.148977f, -1.141770f, -1.134598f, -1.127459f,
98-1.120354f, -1.113282f, -1.106242f, -1.099233f, -1.092255f, // .100
99
100-1.085306f, -1.078388f, -1.071498f, -1.064636f, -1.057802f,
101-1.050996f, -1.044215f, -1.037461f, -1.030733f, -1.024029f,
102-1.017350f, -1.010695f, -1.004064f, -.997456f, -.990871f,
103-.984308f, -.977767f, -.971247f, -.964749f, -.958271f,
104-.951813f, -.945375f, -.938957f, -.932558f, -.926178f,
105-.919816f, -.913472f, -.907146f, -.900838f, -.894547f,
106-.888272f, -.882014f, -.875773f, -.869547f, -.863337f,
107-.857142f, -.850963f, -.844798f, -.838648f, -.832512f,
108-.826390f, -.820282f, -.814187f, -.808106f, -.802038f,
109-.795982f, -.789940f, -.783909f, -.777891f, -.771884f, // .150
110-.765889f, -.759906f, -.753934f, -.747973f, -.742023f,
111-.736084f, -.730155f, -.724237f, -.718328f, -.712429f,
112-.706541f, -.700661f, -.694791f, -.688931f, -.683079f,
113-.677236f, -.671402f, -.665576f, -.659759f, -.653950f,
114-.648149f, -.642356f, -.636570f, -.630793f, -.625022f,
115-.619259f, -.613503f, -.607754f, -.602012f, -.596276f,
116-.590548f, -.584825f, -.579109f, -.573399f, -.567695f,
117-.561997f, -.556305f, -.550618f, -.544937f, -.539262f,
118-.533592f, -.527926f, -.522266f, -.516611f, -.510961f,
119-.505315f, -.499674f, -.494037f, -.488405f, -.482777f, // .200
120
121-.477153f, -.471533f, -.465917f, -.460305f, -.454697f,
122-.449092f, -.443491f, -.437893f, -.432299f, -.426707f,
123-.421119f, -.415534f, -.409951f, -.404372f, -.398795f,
124-.393221f, -.387649f, -.382080f, -.376513f, -.370949f,
125-.365387f, -.359826f, -.354268f, -.348712f, -.343157f,
126-.337604f, -.332053f, -.326503f, -.320955f, -.315408f,
127-.309863f, -.304318f, -.298775f, -.293233f, -.287692f,
128-.282152f, -.276613f, -.271074f, -.265536f, -.259999f,
129-.254462f, -.248926f, -.243389f, -.237854f, -.232318f,
130-.226783f, -.221247f, -.215712f, -.210176f, -.204641f, // .250
131-.199105f, -.193568f, -.188032f, -.182495f, -.176957f,
132-.171419f, -.165880f, -.160341f, -.154800f, -.149259f,
133-.143717f, -.138173f, -.132629f, -.127083f, -.121537f,
134-.115989f, -.110439f, -.104889f, -.099336f, -.093782f,
135-.088227f, -.082670f, -.077111f, -.071550f, -.065987f,
136-.060423f, -.054856f, -.049288f, -.043717f, -.038144f,
137-.032569f, -.026991f, -.021411f, -.015828f, -.010243f,
138-.004656f, .000934f, .006527f, .012123f, .017722f,
139.023323f, .028928f, .034535f, .040146f, .045759f,
140.051376f, .056997f, .062620f, .068247f, .073877f, // .300
141
142.079511f, .085149f, .090790f, .096435f, .102083f,
143.107736f, .113392f, .119052f, .124716f, .130385f,
144.136057f, .141734f, .147414f, .153100f, .158789f,
145.164483f, .170181f, .175884f, .181592f, .187304f,
146.193021f, .198743f, .204469f, .210201f, .215937f,
147.221678f, .227425f, .233177f, .238933f, .244696f,
148.250463f, .256236f, .262014f, .267798f, .273587f,
149.279382f, .285183f, .290989f, .296801f, .302619f,
150.308443f, .314273f, .320109f, .325951f, .331799f,
151.337654f, .343515f, .349382f, .355255f, .361135f, // .350
152.367022f, .372915f, .378815f, .384721f, .390634f,
153.396554f, .402481f, .408415f, .414356f, .420304f,
154.426260f, .432222f, .438192f, .444169f, .450153f,
155.456145f, .462144f, .468151f, .474166f, .480188f,
156.486218f, .492256f, .498302f, .504356f, .510418f,
157.516488f, .522566f, .528653f, .534747f, .540850f,
158.546962f, .553082f, .559210f, .565347f, .571493f,
159.577648f, .583811f, .589983f, .596164f, .602355f,
160.608554f, .614762f, .620980f, .627207f, .633444f,
161.639689f, .645945f, .652210f, .658484f, .664768f, // .400
162
163.671062f, .677366f, .683680f, .690004f, .696338f,
164.702682f, .709036f, .715400f, .721775f, .728160f,
165.734556f, .740963f, .747379f, .753807f, .760246f,
166.766695f, .773155f, .779627f, .786109f, .792603f,
167.799107f, .805624f, .812151f, .818690f, .825241f,
168.831803f, .838377f, .844962f, .851560f, .858170f,
169.864791f, .871425f, .878071f, .884729f, .891399f,
170.898082f, .904778f, .911486f, .918206f, .924940f,
171.931686f, .938446f, .945218f, .952003f, .958802f,
172.965614f, .972439f, .979278f, .986130f, .992996f, // .450
173.999875f, 1.006769f, 1.013676f, 1.020597f, 1.027533f,
1741.034482f, 1.041446f, 1.048424f, 1.055417f, 1.062424f,
1751.069446f, 1.076482f, 1.083534f, 1.090600f, 1.097681f,
1761.104778f, 1.111889f, 1.119016f, 1.126159f, 1.133316f,
1771.140490f, 1.147679f, 1.154884f, 1.162105f, 1.169342f,
1781.176595f, 1.183864f, 1.191149f, 1.198451f, 1.205770f,
1791.213105f, 1.220457f, 1.227826f, 1.235211f, 1.242614f,
1801.250034f, 1.257471f, 1.264926f, 1.272398f, 1.279888f,
1811.287395f, 1.294921f, 1.302464f, 1.310026f, 1.317605f,
1821.325203f, 1.332819f, 1.340454f, 1.348108f, 1.355780f, // .500
183
1841.363472f, 1.371182f, 1.378912f, 1.386660f, 1.394429f,
1851.402216f, 1.410024f, 1.417851f, 1.425698f, 1.433565f,
1861.441453f, 1.449360f, 1.457288f, 1.465237f, 1.473206f,
1871.481196f, 1.489208f, 1.497240f, 1.505293f, 1.513368f,
1881.521465f, 1.529583f, 1.537723f, 1.545885f, 1.554068f,
1891.562275f, 1.570503f, 1.578754f, 1.587028f, 1.595325f,
1901.603644f, 1.611987f, 1.620353f, 1.628743f, 1.637156f,
1911.645593f, 1.654053f, 1.662538f, 1.671047f, 1.679581f,
1921.688139f, 1.696721f, 1.705329f, 1.713961f, 1.722619f,
1931.731303f, 1.740011f, 1.748746f, 1.757506f, 1.766293f, // .550
1941.775106f, 1.783945f, 1.792810f, 1.801703f, 1.810623f,
1951.819569f, 1.828543f, 1.837545f, 1.846574f, 1.855631f,
1961.864717f, 1.873830f, 1.882972f, 1.892143f, 1.901343f,
1971.910572f, 1.919830f, 1.929117f, 1.938434f, 1.947781f,
1981.957158f, 1.966566f, 1.976004f, 1.985473f, 1.994972f,
1992.004503f, 2.014065f, 2.023659f, 2.033285f, 2.042943f,
2002.052633f, 2.062355f, 2.072110f, 2.081899f, 2.091720f,
2012.101575f, 2.111464f, 2.121386f, 2.131343f, 2.141334f,
2022.151360f, 2.161421f, 2.171517f, 2.181648f, 2.191815f,
2032.202018f, 2.212257f, 2.222533f, 2.232845f, 2.243195f, // .600
204
2052.253582f, 2.264006f, 2.274468f, 2.284968f, 2.295507f,
2062.306084f, 2.316701f, 2.327356f, 2.338051f, 2.348786f,
2072.359562f, 2.370377f, 2.381234f, 2.392131f, 2.403070f,
2082.414051f, 2.425073f, 2.436138f, 2.447246f, 2.458397f,
2092.469591f, 2.480828f, 2.492110f, 2.503436f, 2.514807f,
2102.526222f, 2.537684f, 2.549190f, 2.560743f, 2.572343f,
2112.583989f, 2.595682f, 2.607423f, 2.619212f, 2.631050f,
2122.642936f, 2.654871f, 2.666855f, 2.678890f, 2.690975f,
2132.703110f, 2.715297f, 2.727535f, 2.739825f, 2.752168f,
2142.764563f, 2.777012f, 2.789514f, 2.802070f, 2.814681f, // .650
2152.827347f, 2.840069f, 2.852846f, 2.865680f, 2.878570f,
2162.891518f, 2.904524f, 2.917588f, 2.930712f, 2.943894f,
2172.957136f, 2.970439f, 2.983802f, 2.997227f, 3.010714f,
2183.024263f, 3.037875f, 3.051551f, 3.065290f, 3.079095f,
2193.092965f, 3.106900f, 3.120902f, 3.134971f, 3.149107f,
2203.163312f, 3.177585f, 3.191928f, 3.206340f, 3.220824f,
2213.235378f, 3.250005f, 3.264704f, 3.279477f, 3.294323f,
2223.309244f, 3.324240f, 3.339312f, 3.354461f, 3.369687f,
2233.384992f, 3.400375f, 3.415838f, 3.431381f, 3.447005f,
2243.462711f, 3.478500f, 3.494372f, 3.510328f, 3.526370f, // .700
225
2263.542497f, 3.558711f, 3.575012f, 3.591402f, 3.607881f,
2273.624450f, 3.641111f, 3.657863f, 3.674708f, 3.691646f,
2283.708680f, 3.725809f, 3.743034f, 3.760357f, 3.777779f,
2293.795300f, 3.812921f, 3.830645f, 3.848470f, 3.866400f,
2303.884434f, 3.902574f, 3.920821f, 3.939176f, 3.957640f,
2313.976215f, 3.994901f, 4.013699f, 4.032612f, 4.051639f,
2324.070783f, 4.090045f, 4.109425f, 4.128925f, 4.148547f,
2334.168292f, 4.188160f, 4.208154f, 4.228275f, 4.248524f,
2344.268903f, 4.289413f, 4.310056f, 4.330832f, 4.351745f,
2354.372794f, 4.393982f, 4.415310f, 4.436781f, 4.458395f,
2364.480154f, 4.502060f, 4.524114f, 4.546319f, 4.568676f, // .750
2374.591187f, 4.613854f, 4.636678f, 4.659662f, 4.682807f,
2384.706116f, 4.729590f, 4.753231f, 4.777041f, 4.801024f,
2394.825179f, 4.849511f, 4.874020f, 4.898710f, 4.923582f,
2404.948639f, 4.973883f, 4.999316f, 5.024942f, 5.050761f,
2415.076778f, 5.102993f, 5.129411f, 5.156034f, 5.182864f,
2425.209903f, 5.237156f, 5.264625f, 5.292312f, 5.320220f,
2435.348354f, 5.376714f, 5.405306f, 5.434131f, 5.463193f,
2445.492496f, 5.522042f, 5.551836f, 5.581880f, 5.612178f,
2455.642734f, 5.673552f, 5.704634f, 5.735986f, 5.767610f, // .800
246
2475.799512f, 5.831694f, 5.864161f, 5.896918f, 5.929968f,
2485.963316f, 5.996967f, 6.030925f, 6.065194f, 6.099780f,
2496.134687f, 6.169921f, 6.205486f, 6.241387f, 6.277630f,
2506.314220f, 6.351163f, 6.388465f, 6.426130f, 6.464166f,
2516.502578f, 6.541371f, 6.580553f, 6.620130f, 6.660109f,
2526.700495f, 6.741297f, 6.782520f, 6.824173f, 6.866262f,
2536.908795f, 6.951780f, 6.995225f, 7.039137f, 7.083525f,
2547.128398f, 7.173764f, 7.219632f, 7.266011f, 7.312910f,
2557.360339f, 7.408308f, 7.456827f, 7.505905f, 7.555554f,
2567.605785f, 7.656608f, 7.708035f, 7.760077f, 7.812747f, // .850
2577.866057f, 7.920019f, 7.974647f, 8.029953f, 8.085952f,
2588.142657f, 8.200083f, 8.258245f, 8.317158f, 8.376837f,
2598.437300f, 8.498562f, 8.560641f, 8.623554f, 8.687319f,
2608.751955f, 8.817481f, 8.883916f, 8.951282f, 9.019600f,
2619.088889f, 9.159174f, 9.230477f, 9.302822f, 9.376233f,
2629.450735f, 9.526355f, 9.603118f, 9.681054f, 9.760191f,
263 9.840558f, 9.922186f, 10.005107f, 10.089353f, 10.174959f,
26410.261958f, 10.350389f, 10.440287f, 10.531693f, 10.624646f,
26510.719188f, 10.815362f, 10.913214f, 11.012789f, 11.114137f,
26611.217307f, 11.322352f, 11.429325f, 11.538283f, 11.649285f, // .900
267
26811.762390f, 11.877664f, 11.995170f, 12.114979f, 12.237161f,
26912.361791f, 12.488946f, 12.618708f, 12.751161f, 12.886394f,
27013.024498f, 13.165570f, 13.309711f, 13.457026f, 13.607625f,
27113.761625f, 13.919145f, 14.080314f, 14.245263f, 14.414134f,
27214.587072f, 14.764233f, 14.945778f, 15.131877f, 15.322712f,
27315.518470f, 15.719353f, 15.925570f, 16.137345f, 16.354912f,
27416.578520f, 16.808433f, 17.044929f, 17.288305f, 17.538873f,
27517.796967f, 18.062943f, 18.337176f, 18.620068f, 18.912049f,
27619.213574f, 19.525133f, 19.847249f, 20.180480f, 20.525429f,
27720.882738f, 21.253102f, 21.637266f, 22.036036f, 22.450278f, // .950
27822.880933f, 23.329017f, 23.795634f, 24.281981f, 24.789364f,
27925.319207f, 25.873062f, 26.452634f, 27.059789f, 27.696581f, // .960
28028.365274f, 29.068370f, 29.808638f, 30.589157f, 31.413354f,
28132.285060f, 33.208568f, 34.188705f, 35.230920f, 36.341388f, // .970
28237.527131f, 38.796172f, 40.157721f, 41.622399f, 43.202525f,
28344.912465f, 46.769077f, 48.792279f, 51.005773f, 53.437996f, // .980
28456.123356f, 59.103894f, // .982
285
286}; // End of the inverseLandau table
287
288
289double RandLandau::transform (double r) {
290
291 double u = r * TABLE_MULTIPLIER;
292 int index = int(u);
293 double du = u - index;
294
295 // du is scaled such that the we dont have to multiply by TABLE_INTERVAL
296 // when interpolating.
297
298 // Five cases:
299 // A) Between .070 and .800 the function is so smooth, straight
300 // linear interpolation is adequate.
301 // B) Between .007 and .070, and between .800 and .980, quadratic
302 // interpolation is used. This requires the same 4 points as
303 // a cubic spline (thus we need .006 and .981 and .982) but
304 // the quadratic interpolation is accurate enough and quicker.
305 // C) Below .007 an asymptotic expansion for low negative lambda
306 // (involving two logs) is used; there is a pade-style correction
307 // factor.
308 // D) Above .980, a simple pade approximation is made (asymptotic to
309 // 1/(1-r)), but...
310 // E) the coefficients in that pade are different above r=.999.
311
312 if ( index >= 70 && index <= 800 ) { // (A)
313
314 double f0 = inverseLandau [index];
315 double f1 = inverseLandau [index+1];
316 return f0 + du * (f1 - f0);
317
318 } else if ( index >= 7 && index <= 980 ) { // (B)
319
320 double f_1 = inverseLandau [index-1];
321 double f0 = inverseLandau [index];
322 double f1 = inverseLandau [index+1];
323 double f2 = inverseLandau [index+2];
324
325 return f0 + du * (f1 - f0 - .25*(1-du)* (f2 -f1 - f0 + f_1) );
326
327 } else if ( index < 7 ) { // (C)
328
329 const double n0 = 0.99858950;
330 const double n1 = 34.5213058; const double d1 = 34.1760202;
331 const double n2 = 17.0854528; const double d2 = 4.01244582;
332
333 double logr = std::log(r);
334 double x = 1/logr;
335 double x2 = x*x;
336
337 double pade = (n0 + n1*x + n2*x2) / (1.0 + d1*x + d2*x2);
338
339 return ( - std::log ( -.91893853 - logr ) -1 ) * pade;
340
341 } else if ( index <= 999 ) { // (D)
342
343 const double n0 = 1.00060006;
344 const double n1 = 263.991156; const double d1 = 257.368075;
345 const double n2 = 4373.20068; const double d2 = 3414.48018;
346
347 double x = 1-r;
348 double x2 = x*x;
349
350 return (n0 + n1*x + n2*x2) / (x * (1.0 + d1*x + d2*x2));
351
352 } else { // (E)
353
354 const double n0 = 1.00001538;
355 const double n1 = 6075.14119; const double d1 = 6065.11919;
356 const double n2 = 734266.409; const double d2 = 694021.044;
357
358 double x = 1-r;
359 double x2 = x*x;
360
361 return (n0 + n1*x + n2*x2) / (x * (1.0 + d1*x + d2*x2));
362
363 }
364
365} // transform()
366
367std::ostream & RandLandau::put ( std::ostream & os ) const {
368 int pr=os.precision(20);
369 os << " " << name() << "\n";
370 os.precision(pr);
371 return os;
372}
373
374std::istream & RandLandau::get ( std::istream & is ) {
375 std::string inName;
376 is >> inName;
377 if (inName != name()) {
378 is.clear(std::ios::badbit | is.rdstate());
379 std::cerr << "Mismatch when expecting to read state of a "
380 << name() << " distribution\n"
381 << "Name found was " << inName
382 << "\nistream is left in the badbit state\n";
383 return is;
384 }
385 return is;
386}
387
388} // namespace CLHEP
static double shoot()
std::string name() const
Definition RandLandau.cc:28
static double transform(double r)
static void shootArray(const int size, double *vect)
Definition RandLandau.cc:34
virtual ~RandLandau()
Definition RandLandau.cc:31
std::istream & get(std::istream &is)
void fireArray(const int size, double *vect)
Definition RandLandau.cc:48
std::ostream & put(std::ostream &os) const
HepRandomEngine & engine()
Definition RandLandau.cc:29
int(* f2)(int)
void(* f1)()