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