ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RandLandau.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RandLandau.cc
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  int 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