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"
22 #include "CLHEP/Random/RandLandau.h"
23 #include <iostream>
24 #include <cmath> // for std::log()
25 
26 namespace CLHEP {
27 
28 std::string RandLandau::name() const {return "RandLandau";}
29 HepRandomEngine & RandLandau::engine() {return *localEngine;}
30 
32 }
33 
34 void 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 
48 void 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 
62 static const float TABLE_INTERVAL = .001f;
63 static const int TABLE_END = 982;
64 static 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 
76 static const float inverseLandau [TABLE_END+1] = {
77 
78 0.0f, // .000
79 0.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,
174 1.034482f, 1.041446f, 1.048424f, 1.055417f, 1.062424f,
175 1.069446f, 1.076482f, 1.083534f, 1.090600f, 1.097681f,
176 1.104778f, 1.111889f, 1.119016f, 1.126159f, 1.133316f,
177 1.140490f, 1.147679f, 1.154884f, 1.162105f, 1.169342f,
178 1.176595f, 1.183864f, 1.191149f, 1.198451f, 1.205770f,
179 1.213105f, 1.220457f, 1.227826f, 1.235211f, 1.242614f,
180 1.250034f, 1.257471f, 1.264926f, 1.272398f, 1.279888f,
181 1.287395f, 1.294921f, 1.302464f, 1.310026f, 1.317605f,
182 1.325203f, 1.332819f, 1.340454f, 1.348108f, 1.355780f, // .500
183 
184 1.363472f, 1.371182f, 1.378912f, 1.386660f, 1.394429f,
185 1.402216f, 1.410024f, 1.417851f, 1.425698f, 1.433565f,
186 1.441453f, 1.449360f, 1.457288f, 1.465237f, 1.473206f,
187 1.481196f, 1.489208f, 1.497240f, 1.505293f, 1.513368f,
188 1.521465f, 1.529583f, 1.537723f, 1.545885f, 1.554068f,
189 1.562275f, 1.570503f, 1.578754f, 1.587028f, 1.595325f,
190 1.603644f, 1.611987f, 1.620353f, 1.628743f, 1.637156f,
191 1.645593f, 1.654053f, 1.662538f, 1.671047f, 1.679581f,
192 1.688139f, 1.696721f, 1.705329f, 1.713961f, 1.722619f,
193 1.731303f, 1.740011f, 1.748746f, 1.757506f, 1.766293f, // .550
194 1.775106f, 1.783945f, 1.792810f, 1.801703f, 1.810623f,
195 1.819569f, 1.828543f, 1.837545f, 1.846574f, 1.855631f,
196 1.864717f, 1.873830f, 1.882972f, 1.892143f, 1.901343f,
197 1.910572f, 1.919830f, 1.929117f, 1.938434f, 1.947781f,
198 1.957158f, 1.966566f, 1.976004f, 1.985473f, 1.994972f,
199 2.004503f, 2.014065f, 2.023659f, 2.033285f, 2.042943f,
200 2.052633f, 2.062355f, 2.072110f, 2.081899f, 2.091720f,
201 2.101575f, 2.111464f, 2.121386f, 2.131343f, 2.141334f,
202 2.151360f, 2.161421f, 2.171517f, 2.181648f, 2.191815f,
203 2.202018f, 2.212257f, 2.222533f, 2.232845f, 2.243195f, // .600
204 
205 2.253582f, 2.264006f, 2.274468f, 2.284968f, 2.295507f,
206 2.306084f, 2.316701f, 2.327356f, 2.338051f, 2.348786f,
207 2.359562f, 2.370377f, 2.381234f, 2.392131f, 2.403070f,
208 2.414051f, 2.425073f, 2.436138f, 2.447246f, 2.458397f,
209 2.469591f, 2.480828f, 2.492110f, 2.503436f, 2.514807f,
210 2.526222f, 2.537684f, 2.549190f, 2.560743f, 2.572343f,
211 2.583989f, 2.595682f, 2.607423f, 2.619212f, 2.631050f,
212 2.642936f, 2.654871f, 2.666855f, 2.678890f, 2.690975f,
213 2.703110f, 2.715297f, 2.727535f, 2.739825f, 2.752168f,
214 2.764563f, 2.777012f, 2.789514f, 2.802070f, 2.814681f, // .650
215 2.827347f, 2.840069f, 2.852846f, 2.865680f, 2.878570f,
216 2.891518f, 2.904524f, 2.917588f, 2.930712f, 2.943894f,
217 2.957136f, 2.970439f, 2.983802f, 2.997227f, 3.010714f,
218 3.024263f, 3.037875f, 3.051551f, 3.065290f, 3.079095f,
219 3.092965f, 3.106900f, 3.120902f, 3.134971f, 3.149107f,
220 3.163312f, 3.177585f, 3.191928f, 3.206340f, 3.220824f,
221 3.235378f, 3.250005f, 3.264704f, 3.279477f, 3.294323f,
222 3.309244f, 3.324240f, 3.339312f, 3.354461f, 3.369687f,
223 3.384992f, 3.400375f, 3.415838f, 3.431381f, 3.447005f,
224 3.462711f, 3.478500f, 3.494372f, 3.510328f, 3.526370f, // .700
225 
226 3.542497f, 3.558711f, 3.575012f, 3.591402f, 3.607881f,
227 3.624450f, 3.641111f, 3.657863f, 3.674708f, 3.691646f,
228 3.708680f, 3.725809f, 3.743034f, 3.760357f, 3.777779f,
229 3.795300f, 3.812921f, 3.830645f, 3.848470f, 3.866400f,
230 3.884434f, 3.902574f, 3.920821f, 3.939176f, 3.957640f,
231 3.976215f, 3.994901f, 4.013699f, 4.032612f, 4.051639f,
232 4.070783f, 4.090045f, 4.109425f, 4.128925f, 4.148547f,
233 4.168292f, 4.188160f, 4.208154f, 4.228275f, 4.248524f,
234 4.268903f, 4.289413f, 4.310056f, 4.330832f, 4.351745f,
235 4.372794f, 4.393982f, 4.415310f, 4.436781f, 4.458395f,
236 4.480154f, 4.502060f, 4.524114f, 4.546319f, 4.568676f, // .750
237 4.591187f, 4.613854f, 4.636678f, 4.659662f, 4.682807f,
238 4.706116f, 4.729590f, 4.753231f, 4.777041f, 4.801024f,
239 4.825179f, 4.849511f, 4.874020f, 4.898710f, 4.923582f,
240 4.948639f, 4.973883f, 4.999316f, 5.024942f, 5.050761f,
241 5.076778f, 5.102993f, 5.129411f, 5.156034f, 5.182864f,
242 5.209903f, 5.237156f, 5.264625f, 5.292312f, 5.320220f,
243 5.348354f, 5.376714f, 5.405306f, 5.434131f, 5.463193f,
244 5.492496f, 5.522042f, 5.551836f, 5.581880f, 5.612178f,
245 5.642734f, 5.673552f, 5.704634f, 5.735986f, 5.767610f, // .800
246 
247 5.799512f, 5.831694f, 5.864161f, 5.896918f, 5.929968f,
248 5.963316f, 5.996967f, 6.030925f, 6.065194f, 6.099780f,
249 6.134687f, 6.169921f, 6.205486f, 6.241387f, 6.277630f,
250 6.314220f, 6.351163f, 6.388465f, 6.426130f, 6.464166f,
251 6.502578f, 6.541371f, 6.580553f, 6.620130f, 6.660109f,
252 6.700495f, 6.741297f, 6.782520f, 6.824173f, 6.866262f,
253 6.908795f, 6.951780f, 6.995225f, 7.039137f, 7.083525f,
254 7.128398f, 7.173764f, 7.219632f, 7.266011f, 7.312910f,
255 7.360339f, 7.408308f, 7.456827f, 7.505905f, 7.555554f,
256 7.605785f, 7.656608f, 7.708035f, 7.760077f, 7.812747f, // .850
257 7.866057f, 7.920019f, 7.974647f, 8.029953f, 8.085952f,
258 8.142657f, 8.200083f, 8.258245f, 8.317158f, 8.376837f,
259 8.437300f, 8.498562f, 8.560641f, 8.623554f, 8.687319f,
260 8.751955f, 8.817481f, 8.883916f, 8.951282f, 9.019600f,
261 9.088889f, 9.159174f, 9.230477f, 9.302822f, 9.376233f,
262 9.450735f, 9.526355f, 9.603118f, 9.681054f, 9.760191f,
263  9.840558f, 9.922186f, 10.005107f, 10.089353f, 10.174959f,
264 10.261958f, 10.350389f, 10.440287f, 10.531693f, 10.624646f,
265 10.719188f, 10.815362f, 10.913214f, 11.012789f, 11.114137f,
266 11.217307f, 11.322352f, 11.429325f, 11.538283f, 11.649285f, // .900
267 
268 11.762390f, 11.877664f, 11.995170f, 12.114979f, 12.237161f,
269 12.361791f, 12.488946f, 12.618708f, 12.751161f, 12.886394f,
270 13.024498f, 13.165570f, 13.309711f, 13.457026f, 13.607625f,
271 13.761625f, 13.919145f, 14.080314f, 14.245263f, 14.414134f,
272 14.587072f, 14.764233f, 14.945778f, 15.131877f, 15.322712f,
273 15.518470f, 15.719353f, 15.925570f, 16.137345f, 16.354912f,
274 16.578520f, 16.808433f, 17.044929f, 17.288305f, 17.538873f,
275 17.796967f, 18.062943f, 18.337176f, 18.620068f, 18.912049f,
276 19.213574f, 19.525133f, 19.847249f, 20.180480f, 20.525429f,
277 20.882738f, 21.253102f, 21.637266f, 22.036036f, 22.450278f, // .950
278 22.880933f, 23.329017f, 23.795634f, 24.281981f, 24.789364f,
279 25.319207f, 25.873062f, 26.452634f, 27.059789f, 27.696581f, // .960
280 28.365274f, 29.068370f, 29.808638f, 30.589157f, 31.413354f,
281 32.285060f, 33.208568f, 34.188705f, 35.230920f, 36.341388f, // .970
282 37.527131f, 38.796172f, 40.157721f, 41.622399f, 43.202525f,
283 44.912465f, 46.769077f, 48.792279f, 51.005773f, 53.437996f, // .980
284 56.123356f, 59.103894f, // .982
285 
286 }; // End of the inverseLandau table
287 
288 
289 double 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 
367 std::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 
374 std::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
HepRandomEngine & engine()
Definition: RandLandau.cc:29
static void shootArray(const int size, double *vect)
Definition: RandLandau.cc:34
void(* f1)()
void fireArray(const int size, double *vect)
Definition: RandLandau.cc:48
std::string name() const
Definition: RandLandau.cc:28
int(* f2)(int)
std::ostream & put(std::ostream &os) const
Definition: RandLandau.cc:367
static double shoot()
static double transform(double r)
Definition: RandLandau.cc:289
virtual ~RandLandau()
Definition: RandLandau.cc:31
std::istream & get(std::istream &is)
Definition: RandLandau.cc:374