| Geant4 Cross Reference |
1 // -*- C++ -*-
2 //
3 // -----------------------------------------------------------------------
4 // HEP Random
5 // --- RandLandau ---
6 // class implementation file
7 // -----------------------------------------------------------------------
8
9 // =======================================================================
10 // M Fischler - Created 1/6/2000.
11 //
12 // The key transform() method uses the algorithm in CERNLIB.
13 // This is because I trust that RANLAN routine more than
14 // I trust the Bukin-Grozina inverseLandau, which is not
15 // claimed to be better than 1% accurate.
16 //
17 // M Fischler - put and get to/from streams 12/13/04
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 "RandLandau";}
27 HepRandomEngine & RandLandau::engine() {return *localEngine;}
28
29 RandLandau::~RandLandau() {
30 }
31
32 void RandLandau::shootArray( const int size, double* vect )
33
34 {
35 for( double* v = vect; v != vect + size; ++v )
36 *v = shoot();
37 }
38
39 void RandLandau::shootArray( HepRandomEngine* anEngine,
40 const int size, double* vect )
41 {
42 for( double* v = vect; v != vect + size; ++v )
43 *v = shoot(anEngine);
44 }
45
46 void RandLandau::fireArray( const int size, double* vect)
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 = .060 to .982
54 //
55
56 // Since all these are this is static to this compilation unit only, the
57 // info is establised a priori and not at each invocation.
58
59 static const float TABLE_INTERVAL = .001f;
60 static const int TABLE_END = 982;
61 static const float TABLE_MULTIPLIER = 1.0f/TABLE_INTERVAL;
62
63 // Here comes the big (4K bytes) table ---
64 //
65 // inverseLandau[ n ] = the inverse cdf at r = n*TABLE_INTERVAL = n/1000.
66 //
67 // Credit CERNLIB for these computations.
68 //
69 // This data is float because the algortihm does not benefit from further
70 // data accuracy. The numbers below .006 or above .982 are moot since
71 // non-table-based methods are used below r=.007 and above .980.
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, // .001 - .005
77 -2.244733f, -2.204365f, -2.168163f, -2.135219f, -2.104898f, // .006 - .010
78 -2.076740f, -2.050397f, -2.025605f, -2.002150f, -1.979866f,
79 -1.958612f, -1.938275f, -1.918760f, -1.899984f, -1.881879f, // .020
80 -1.864385f, -1.847451f, -1.831030f, -1.815083f, -1.799574f,
81 -1.784473f, -1.769751f, -1.755383f, -1.741346f, -1.727620f, // .030
82 -1.714187f, -1.701029f, -1.688130f, -1.675477f, -1.663057f,
83 -1.650858f, -1.638868f, -1.627078f, -1.615477f, -1.604058f, // .040
84 -1.592811f, -1.581729f, -1.570806f, -1.560034f, -1.549407f,
85 -1.538919f, -1.528565f, -1.518339f, -1.508237f, -1.498254f, // .050
86 -1.488386f, -1.478628f, -1.468976f, -1.459428f, -1.449979f,
87 -1.440626f, -1.431365f, -1.422195f, -1.413111f, -1.404112f, // .060
88 -1.395194f, -1.386356f, -1.377594f, -1.368906f, -1.360291f,
89 -1.351746f, -1.343269f, -1.334859f, -1.326512f, -1.318229f, // .070
90 -1.310006f, -1.301843f, -1.293737f, -1.285688f, -1.277693f,
91 -1.269752f, -1.261863f, -1.254024f, -1.246235f, -1.238494f, // .080
92 -1.230800f, -1.223153f, -1.215550f, -1.207990f, -1.200474f,
93 -1.192999f, -1.185566f, -1.178172f, -1.170817f, -1.163500f, // .090
94 -1.156220f, -1.148977f, -1.141770f, -1.134598f, -1.127459f,
95 -1.120354f, -1.113282f, -1.106242f, -1.099233f, -1.092255f, // .100
96
97 -1.085306f, -1.078388f, -1.071498f, -1.064636f, -1.057802f,
98 -1.050996f, -1.044215f, -1.037461f, -1.030733f, -1.024029f,
99 -1.017350f, -1.010695f, -1.004064f, -.997456f, -.990871f,
100 -.984308f, -.977767f, -.971247f, -.964749f, -.958271f,
101 -.951813f, -.945375f, -.938957f, -.932558f, -.926178f,
102 -.919816f, -.913472f, -.907146f, -.900838f, -.894547f,
103 -.888272f, -.882014f, -.875773f, -.869547f, -.863337f,
104 -.857142f, -.850963f, -.844798f, -.838648f, -.832512f,
105 -.826390f, -.820282f, -.814187f, -.808106f, -.802038f,
106 -.795982f, -.789940f, -.783909f, -.777891f, -.771884f, // .150
107 -.765889f, -.759906f, -.753934f, -.747973f, -.742023f,
108 -.736084f, -.730155f, -.724237f, -.718328f, -.712429f,
109 -.706541f, -.700661f, -.694791f, -.688931f, -.683079f,
110 -.677236f, -.671402f, -.665576f, -.659759f, -.653950f,
111 -.648149f, -.642356f, -.636570f, -.630793f, -.625022f,
112 -.619259f, -.613503f, -.607754f, -.602012f, -.596276f,
113 -.590548f, -.584825f, -.579109f, -.573399f, -.567695f,
114 -.561997f, -.556305f, -.550618f, -.544937f, -.539262f,
115 -.533592f, -.527926f, -.522266f, -.516611f, -.510961f,
116 -.505315f, -.499674f, -.494037f, -.488405f, -.482777f, // .200
117
118 -.477153f, -.471533f, -.465917f, -.460305f, -.454697f,
119 -.449092f, -.443491f, -.437893f, -.432299f, -.426707f,
120 -.421119f, -.415534f, -.409951f, -.404372f, -.398795f,
121 -.393221f, -.387649f, -.382080f, -.376513f, -.370949f,
122 -.365387f, -.359826f, -.354268f, -.348712f, -.343157f,
123 -.337604f, -.332053f, -.326503f, -.320955f, -.315408f,
124 -.309863f, -.304318f, -.298775f, -.293233f, -.287692f,
125 -.282152f, -.276613f, -.271074f, -.265536f, -.259999f,
126 -.254462f, -.248926f, -.243389f, -.237854f, -.232318f,
127 -.226783f, -.221247f, -.215712f, -.210176f, -.204641f, // .250
128 -.199105f, -.193568f, -.188032f, -.182495f, -.176957f,
129 -.171419f, -.165880f, -.160341f, -.154800f, -.149259f,
130 -.143717f, -.138173f, -.132629f, -.127083f, -.121537f,
131 -.115989f, -.110439f, -.104889f, -.099336f, -.093782f,
132 -.088227f, -.082670f, -.077111f, -.071550f, -.065987f,
133 -.060423f, -.054856f, -.049288f, -.043717f, -.038144f,
134 -.032569f, -.026991f, -.021411f, -.015828f, -.010243f,
135 -.004656f, .000934f, .006527f, .012123f, .017722f,
136 .023323f, .028928f, .034535f, .040146f, .045759f,
137 .051376f, .056997f, .062620f, .068247f, .073877f, // .300
138
139 .079511f, .085149f, .090790f, .096435f, .102083f,
140 .107736f, .113392f, .119052f, .124716f, .130385f,
141 .136057f, .141734f, .147414f, .153100f, .158789f,
142 .164483f, .170181f, .175884f, .181592f, .187304f,
143 .193021f, .198743f, .204469f, .210201f, .215937f,
144 .221678f, .227425f, .233177f, .238933f, .244696f,
145 .250463f, .256236f, .262014f, .267798f, .273587f,
146 .279382f, .285183f, .290989f, .296801f, .302619f,
147 .308443f, .314273f, .320109f, .325951f, .331799f,
148 .337654f, .343515f, .349382f, .355255f, .361135f, // .350
149 .367022f, .372915f, .378815f, .384721f, .390634f,
150 .396554f, .402481f, .408415f, .414356f, .420304f,
151 .426260f, .432222f, .438192f, .444169f, .450153f,
152 .456145f, .462144f, .468151f, .474166f, .480188f,
153 .486218f, .492256f, .498302f, .504356f, .510418f,
154 .516488f, .522566f, .528653f, .534747f, .540850f,
155 .546962f, .553082f, .559210f, .565347f, .571493f,
156 .577648f, .583811f, .589983f, .596164f, .602355f,
157 .608554f, .614762f, .620980f, .627207f, .633444f,
158 .639689f, .645945f, .652210f, .658484f, .664768f, // .400
159
160 .671062f, .677366f, .683680f, .690004f, .696338f,
161 .702682f, .709036f, .715400f, .721775f, .728160f,
162 .734556f, .740963f, .747379f, .753807f, .760246f,
163 .766695f, .773155f, .779627f, .786109f, .792603f,
164 .799107f, .805624f, .812151f, .818690f, .825241f,
165 .831803f, .838377f, .844962f, .851560f, .858170f,
166 .864791f, .871425f, .878071f, .884729f, .891399f,
167 .898082f, .904778f, .911486f, .918206f, .924940f,
168 .931686f, .938446f, .945218f, .952003f, .958802f,
169 .965614f, .972439f, .979278f, .986130f, .992996f, // .450
170 .999875f, 1.006769f, 1.013676f, 1.020597f, 1.027533f,
171 1.034482f, 1.041446f, 1.048424f, 1.055417f, 1.062424f,
172 1.069446f, 1.076482f, 1.083534f, 1.090600f, 1.097681f,
173 1.104778f, 1.111889f, 1.119016f, 1.126159f, 1.133316f,
174 1.140490f, 1.147679f, 1.154884f, 1.162105f, 1.169342f,
175 1.176595f, 1.183864f, 1.191149f, 1.198451f, 1.205770f,
176 1.213105f, 1.220457f, 1.227826f, 1.235211f, 1.242614f,
177 1.250034f, 1.257471f, 1.264926f, 1.272398f, 1.279888f,
178 1.287395f, 1.294921f, 1.302464f, 1.310026f, 1.317605f,
179 1.325203f, 1.332819f, 1.340454f, 1.348108f, 1.355780f, // .500
180
181 1.363472f, 1.371182f, 1.378912f, 1.386660f, 1.394429f,
182 1.402216f, 1.410024f, 1.417851f, 1.425698f, 1.433565f,
183 1.441453f, 1.449360f, 1.457288f, 1.465237f, 1.473206f,
184 1.481196f, 1.489208f, 1.497240f, 1.505293f, 1.513368f,
185 1.521465f, 1.529583f, 1.537723f, 1.545885f, 1.554068f,
186 1.562275f, 1.570503f, 1.578754f, 1.587028f, 1.595325f,
187 1.603644f, 1.611987f, 1.620353f, 1.628743f, 1.637156f,
188 1.645593f, 1.654053f, 1.662538f, 1.671047f, 1.679581f,
189 1.688139f, 1.696721f, 1.705329f, 1.713961f, 1.722619f,
190 1.731303f, 1.740011f, 1.748746f, 1.757506f, 1.766293f, // .550
191 1.775106f, 1.783945f, 1.792810f, 1.801703f, 1.810623f,
192 1.819569f, 1.828543f, 1.837545f, 1.846574f, 1.855631f,
193 1.864717f, 1.873830f, 1.882972f, 1.892143f, 1.901343f,
194 1.910572f, 1.919830f, 1.929117f, 1.938434f, 1.947781f,
195 1.957158f, 1.966566f, 1.976004f, 1.985473f, 1.994972f,
196 2.004503f, 2.014065f, 2.023659f, 2.033285f, 2.042943f,
197 2.052633f, 2.062355f, 2.072110f, 2.081899f, 2.091720f,
198 2.101575f, 2.111464f, 2.121386f, 2.131343f, 2.141334f,
199 2.151360f, 2.161421f, 2.171517f, 2.181648f, 2.191815f,
200 2.202018f, 2.212257f, 2.222533f, 2.232845f, 2.243195f, // .600
201
202 2.253582f, 2.264006f, 2.274468f, 2.284968f, 2.295507f,
203 2.306084f, 2.316701f, 2.327356f, 2.338051f, 2.348786f,
204 2.359562f, 2.370377f, 2.381234f, 2.392131f, 2.403070f,
205 2.414051f, 2.425073f, 2.436138f, 2.447246f, 2.458397f,
206 2.469591f, 2.480828f, 2.492110f, 2.503436f, 2.514807f,
207 2.526222f, 2.537684f, 2.549190f, 2.560743f, 2.572343f,
208 2.583989f, 2.595682f, 2.607423f, 2.619212f, 2.631050f,
209 2.642936f, 2.654871f, 2.666855f, 2.678890f, 2.690975f,
210 2.703110f, 2.715297f, 2.727535f, 2.739825f, 2.752168f,
211 2.764563f, 2.777012f, 2.789514f, 2.802070f, 2.814681f, // .650
212 2.827347f, 2.840069f, 2.852846f, 2.865680f, 2.878570f,
213 2.891518f, 2.904524f, 2.917588f, 2.930712f, 2.943894f,
214 2.957136f, 2.970439f, 2.983802f, 2.997227f, 3.010714f,
215 3.024263f, 3.037875f, 3.051551f, 3.065290f, 3.079095f,
216 3.092965f, 3.106900f, 3.120902f, 3.134971f, 3.149107f,
217 3.163312f, 3.177585f, 3.191928f, 3.206340f, 3.220824f,
218 3.235378f, 3.250005f, 3.264704f, 3.279477f, 3.294323f,
219 3.309244f, 3.324240f, 3.339312f, 3.354461f, 3.369687f,
220 3.384992f, 3.400375f, 3.415838f, 3.431381f, 3.447005f,
221 3.462711f, 3.478500f, 3.494372f, 3.510328f, 3.526370f, // .700
222
223 3.542497f, 3.558711f, 3.575012f, 3.591402f, 3.607881f,
224 3.624450f, 3.641111f, 3.657863f, 3.674708f, 3.691646f,
225 3.708680f, 3.725809f, 3.743034f, 3.760357f, 3.777779f,
226 3.795300f, 3.812921f, 3.830645f, 3.848470f, 3.866400f,
227 3.884434f, 3.902574f, 3.920821f, 3.939176f, 3.957640f,
228 3.976215f, 3.994901f, 4.013699f, 4.032612f, 4.051639f,
229 4.070783f, 4.090045f, 4.109425f, 4.128925f, 4.148547f,
230 4.168292f, 4.188160f, 4.208154f, 4.228275f, 4.248524f,
231 4.268903f, 4.289413f, 4.310056f, 4.330832f, 4.351745f,
232 4.372794f, 4.393982f, 4.415310f, 4.436781f, 4.458395f,
233 4.480154f, 4.502060f, 4.524114f, 4.546319f, 4.568676f, // .750
234 4.591187f, 4.613854f, 4.636678f, 4.659662f, 4.682807f,
235 4.706116f, 4.729590f, 4.753231f, 4.777041f, 4.801024f,
236 4.825179f, 4.849511f, 4.874020f, 4.898710f, 4.923582f,
237 4.948639f, 4.973883f, 4.999316f, 5.024942f, 5.050761f,
238 5.076778f, 5.102993f, 5.129411f, 5.156034f, 5.182864f,
239 5.209903f, 5.237156f, 5.264625f, 5.292312f, 5.320220f,
240 5.348354f, 5.376714f, 5.405306f, 5.434131f, 5.463193f,
241 5.492496f, 5.522042f, 5.551836f, 5.581880f, 5.612178f,
242 5.642734f, 5.673552f, 5.704634f, 5.735986f, 5.767610f, // .800
243
244 5.799512f, 5.831694f, 5.864161f, 5.896918f, 5.929968f,
245 5.963316f, 5.996967f, 6.030925f, 6.065194f, 6.099780f,
246 6.134687f, 6.169921f, 6.205486f, 6.241387f, 6.277630f,
247 6.314220f, 6.351163f, 6.388465f, 6.426130f, 6.464166f,
248 6.502578f, 6.541371f, 6.580553f, 6.620130f, 6.660109f,
249 6.700495f, 6.741297f, 6.782520f, 6.824173f, 6.866262f,
250 6.908795f, 6.951780f, 6.995225f, 7.039137f, 7.083525f,
251 7.128398f, 7.173764f, 7.219632f, 7.266011f, 7.312910f,
252 7.360339f, 7.408308f, 7.456827f, 7.505905f, 7.555554f,
253 7.605785f, 7.656608f, 7.708035f, 7.760077f, 7.812747f, // .850
254 7.866057f, 7.920019f, 7.974647f, 8.029953f, 8.085952f,
255 8.142657f, 8.200083f, 8.258245f, 8.317158f, 8.376837f,
256 8.437300f, 8.498562f, 8.560641f, 8.623554f, 8.687319f,
257 8.751955f, 8.817481f, 8.883916f, 8.951282f, 9.019600f,
258 9.088889f, 9.159174f, 9.230477f, 9.302822f, 9.376233f,
259 9.450735f, 9.526355f, 9.603118f, 9.681054f, 9.760191f,
260 9.840558f, 9.922186f, 10.005107f, 10.089353f, 10.174959f,
261 10.261958f, 10.350389f, 10.440287f, 10.531693f, 10.624646f,
262 10.719188f, 10.815362f, 10.913214f, 11.012789f, 11.114137f,
263 11.217307f, 11.322352f, 11.429325f, 11.538283f, 11.649285f, // .900
264
265 11.762390f, 11.877664f, 11.995170f, 12.114979f, 12.237161f,
266 12.361791f, 12.488946f, 12.618708f, 12.751161f, 12.886394f,
267 13.024498f, 13.165570f, 13.309711f, 13.457026f, 13.607625f,
268 13.761625f, 13.919145f, 14.080314f, 14.245263f, 14.414134f,
269 14.587072f, 14.764233f, 14.945778f, 15.131877f, 15.322712f,
270 15.518470f, 15.719353f, 15.925570f, 16.137345f, 16.354912f,
271 16.578520f, 16.808433f, 17.044929f, 17.288305f, 17.538873f,
272 17.796967f, 18.062943f, 18.337176f, 18.620068f, 18.912049f,
273 19.213574f, 19.525133f, 19.847249f, 20.180480f, 20.525429f,
274 20.882738f, 21.253102f, 21.637266f, 22.036036f, 22.450278f, // .950
275 22.880933f, 23.329017f, 23.795634f, 24.281981f, 24.789364f,
276 25.319207f, 25.873062f, 26.452634f, 27.059789f, 27.696581f, // .960
277 28.365274f, 29.068370f, 29.808638f, 30.589157f, 31.413354f,
278 32.285060f, 33.208568f, 34.188705f, 35.230920f, 36.341388f, // .970
279 37.527131f, 38.796172f, 40.157721f, 41.622399f, 43.202525f,
280 44.912465f, 46.769077f, 48.792279f, 51.005773f, 53.437996f, // .980
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 to multiply by TABLE_INTERVAL
292 // when interpolating.
293
294 // Five cases:
295 // A) Between .070 and .800 the function is so smooth, straight
296 // linear interpolation is adequate.
297 // B) Between .007 and .070, and between .800 and .980, quadratic
298 // interpolation is used. This requires the same 4 points as
299 // a cubic spline (thus we need .006 and .981 and .982) but
300 // the quadratic interpolation is accurate enough and quicker.
301 // C) Below .007 an asymptotic expansion for low negative lambda
302 // (involving two logs) is used; there is a pade-style correction
303 // factor.
304 // D) Above .980, a simple pade approximation is made (asymptotic to
305 // 1/(1-r)), but...
306 // E) the coefficients in that pade are different above r=.999.
307
308 if ( index >= 70 && index <= 800 ) { // (A)
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 ) { // (B)
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)* (f2 -f1 - f0 + f_1) );
322
323 } else if ( index < 7 ) { // (C)
324
325 const double n0 = 0.99858950;
326 const double n1 = 34.5213058; const double d1 = 34.1760202;
327 const double n2 = 17.0854528; const double d2 = 4.01244582;
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 + d1*x + d2*x2);
334
335 return ( - std::log ( -.91893853 - logr ) -1 ) * pade;
336
337 } else if ( index <= 999 ) { // (D)
338
339 const double n0 = 1.00060006;
340 const double n1 = 263.991156; const double d1 = 257.368075;
341 const double n2 = 4373.20068; const double d2 = 3414.48018;
342
343 double x = 1-r;
344 double x2 = x*x;
345
346 return (n0 + n1*x + n2*x2) / (x * (1.0 + d1*x + d2*x2));
347
348 } else { // (E)
349
350 const double n0 = 1.00001538;
351 const double n1 = 6075.14119; const double d1 = 6065.11919;
352 const double n2 = 734266.409; const double d2 = 694021.044;
353
354 double x = 1-r;
355 double x2 = x*x;
356
357 return (n0 + n1*x + n2*x2) / (x * (1.0 + d1*x + d2*x2));
358
359 }
360
361 } // transform()
362
363 std::ostream & RandLandau::put ( std::ostream & os ) const {
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 & is ) {
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 read state of a "
376 << name() << " distribution\n"
377 << "Name found was " << inName
378 << "\nistream is left in the badbit state\n";
379 return is;
380 }
381 return is;
382 }
383
384 } // namespace CLHEP
385