00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "CLHEP/Random/RandLandau.h"
00022 #include <iostream>
00023 #include <cmath>
00024
00025 namespace CLHEP {
00026
00027 std::string RandLandau::name() const {return "RandLandau";}
00028 HepRandomEngine & RandLandau::engine() {return *localEngine;}
00029
00030 RandLandau::~RandLandau() {
00031 }
00032
00033 void RandLandau::shootArray( const int size, double* vect )
00034
00035 {
00036 for( double* v = vect; v != vect + size; ++v )
00037 *v = shoot();
00038 }
00039
00040 void RandLandau::shootArray( HepRandomEngine* anEngine,
00041 const int size, double* vect )
00042 {
00043 for( double* v = vect; v != vect + size; ++v )
00044 *v = shoot(anEngine);
00045 }
00046
00047 void RandLandau::fireArray( const int size, double* vect)
00048 {
00049 for( double* v = vect; v != vect + size; ++v )
00050 *v = fire();
00051 }
00052
00053
00054
00055
00056
00057
00058
00059
00060 static const float TABLE_INTERVAL = .001f;
00061 static const int TABLE_END = 982;
00062 static const float TABLE_MULTIPLIER = 1.0f/TABLE_INTERVAL;
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074 static const float inverseLandau [TABLE_END+1] = {
00075
00076 0.0f,
00077 0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
00078 -2.244733f, -2.204365f, -2.168163f, -2.135219f, -2.104898f,
00079 -2.076740f, -2.050397f, -2.025605f, -2.002150f, -1.979866f,
00080 -1.958612f, -1.938275f, -1.918760f, -1.899984f, -1.881879f,
00081 -1.864385f, -1.847451f, -1.831030f, -1.815083f, -1.799574f,
00082 -1.784473f, -1.769751f, -1.755383f, -1.741346f, -1.727620f,
00083 -1.714187f, -1.701029f, -1.688130f, -1.675477f, -1.663057f,
00084 -1.650858f, -1.638868f, -1.627078f, -1.615477f, -1.604058f,
00085 -1.592811f, -1.581729f, -1.570806f, -1.560034f, -1.549407f,
00086 -1.538919f, -1.528565f, -1.518339f, -1.508237f, -1.498254f,
00087 -1.488386f, -1.478628f, -1.468976f, -1.459428f, -1.449979f,
00088 -1.440626f, -1.431365f, -1.422195f, -1.413111f, -1.404112f,
00089 -1.395194f, -1.386356f, -1.377594f, -1.368906f, -1.360291f,
00090 -1.351746f, -1.343269f, -1.334859f, -1.326512f, -1.318229f,
00091 -1.310006f, -1.301843f, -1.293737f, -1.285688f, -1.277693f,
00092 -1.269752f, -1.261863f, -1.254024f, -1.246235f, -1.238494f,
00093 -1.230800f, -1.223153f, -1.215550f, -1.207990f, -1.200474f,
00094 -1.192999f, -1.185566f, -1.178172f, -1.170817f, -1.163500f,
00095 -1.156220f, -1.148977f, -1.141770f, -1.134598f, -1.127459f,
00096 -1.120354f, -1.113282f, -1.106242f, -1.099233f, -1.092255f,
00097
00098 -1.085306f, -1.078388f, -1.071498f, -1.064636f, -1.057802f,
00099 -1.050996f, -1.044215f, -1.037461f, -1.030733f, -1.024029f,
00100 -1.017350f, -1.010695f, -1.004064f, -.997456f, -.990871f,
00101 -.984308f, -.977767f, -.971247f, -.964749f, -.958271f,
00102 -.951813f, -.945375f, -.938957f, -.932558f, -.926178f,
00103 -.919816f, -.913472f, -.907146f, -.900838f, -.894547f,
00104 -.888272f, -.882014f, -.875773f, -.869547f, -.863337f,
00105 -.857142f, -.850963f, -.844798f, -.838648f, -.832512f,
00106 -.826390f, -.820282f, -.814187f, -.808106f, -.802038f,
00107 -.795982f, -.789940f, -.783909f, -.777891f, -.771884f,
00108 -.765889f, -.759906f, -.753934f, -.747973f, -.742023f,
00109 -.736084f, -.730155f, -.724237f, -.718328f, -.712429f,
00110 -.706541f, -.700661f, -.694791f, -.688931f, -.683079f,
00111 -.677236f, -.671402f, -.665576f, -.659759f, -.653950f,
00112 -.648149f, -.642356f, -.636570f, -.630793f, -.625022f,
00113 -.619259f, -.613503f, -.607754f, -.602012f, -.596276f,
00114 -.590548f, -.584825f, -.579109f, -.573399f, -.567695f,
00115 -.561997f, -.556305f, -.550618f, -.544937f, -.539262f,
00116 -.533592f, -.527926f, -.522266f, -.516611f, -.510961f,
00117 -.505315f, -.499674f, -.494037f, -.488405f, -.482777f,
00118
00119 -.477153f, -.471533f, -.465917f, -.460305f, -.454697f,
00120 -.449092f, -.443491f, -.437893f, -.432299f, -.426707f,
00121 -.421119f, -.415534f, -.409951f, -.404372f, -.398795f,
00122 -.393221f, -.387649f, -.382080f, -.376513f, -.370949f,
00123 -.365387f, -.359826f, -.354268f, -.348712f, -.343157f,
00124 -.337604f, -.332053f, -.326503f, -.320955f, -.315408f,
00125 -.309863f, -.304318f, -.298775f, -.293233f, -.287692f,
00126 -.282152f, -.276613f, -.271074f, -.265536f, -.259999f,
00127 -.254462f, -.248926f, -.243389f, -.237854f, -.232318f,
00128 -.226783f, -.221247f, -.215712f, -.210176f, -.204641f,
00129 -.199105f, -.193568f, -.188032f, -.182495f, -.176957f,
00130 -.171419f, -.165880f, -.160341f, -.154800f, -.149259f,
00131 -.143717f, -.138173f, -.132629f, -.127083f, -.121537f,
00132 -.115989f, -.110439f, -.104889f, -.099336f, -.093782f,
00133 -.088227f, -.082670f, -.077111f, -.071550f, -.065987f,
00134 -.060423f, -.054856f, -.049288f, -.043717f, -.038144f,
00135 -.032569f, -.026991f, -.021411f, -.015828f, -.010243f,
00136 -.004656f, .000934f, .006527f, .012123f, .017722f,
00137 .023323f, .028928f, .034535f, .040146f, .045759f,
00138 .051376f, .056997f, .062620f, .068247f, .073877f,
00139
00140 .079511f, .085149f, .090790f, .096435f, .102083f,
00141 .107736f, .113392f, .119052f, .124716f, .130385f,
00142 .136057f, .141734f, .147414f, .153100f, .158789f,
00143 .164483f, .170181f, .175884f, .181592f, .187304f,
00144 .193021f, .198743f, .204469f, .210201f, .215937f,
00145 .221678f, .227425f, .233177f, .238933f, .244696f,
00146 .250463f, .256236f, .262014f, .267798f, .273587f,
00147 .279382f, .285183f, .290989f, .296801f, .302619f,
00148 .308443f, .314273f, .320109f, .325951f, .331799f,
00149 .337654f, .343515f, .349382f, .355255f, .361135f,
00150 .367022f, .372915f, .378815f, .384721f, .390634f,
00151 .396554f, .402481f, .408415f, .414356f, .420304f,
00152 .426260f, .432222f, .438192f, .444169f, .450153f,
00153 .456145f, .462144f, .468151f, .474166f, .480188f,
00154 .486218f, .492256f, .498302f, .504356f, .510418f,
00155 .516488f, .522566f, .528653f, .534747f, .540850f,
00156 .546962f, .553082f, .559210f, .565347f, .571493f,
00157 .577648f, .583811f, .589983f, .596164f, .602355f,
00158 .608554f, .614762f, .620980f, .627207f, .633444f,
00159 .639689f, .645945f, .652210f, .658484f, .664768f,
00160
00161 .671062f, .677366f, .683680f, .690004f, .696338f,
00162 .702682f, .709036f, .715400f, .721775f, .728160f,
00163 .734556f, .740963f, .747379f, .753807f, .760246f,
00164 .766695f, .773155f, .779627f, .786109f, .792603f,
00165 .799107f, .805624f, .812151f, .818690f, .825241f,
00166 .831803f, .838377f, .844962f, .851560f, .858170f,
00167 .864791f, .871425f, .878071f, .884729f, .891399f,
00168 .898082f, .904778f, .911486f, .918206f, .924940f,
00169 .931686f, .938446f, .945218f, .952003f, .958802f,
00170 .965614f, .972439f, .979278f, .986130f, .992996f,
00171 .999875f, 1.006769f, 1.013676f, 1.020597f, 1.027533f,
00172 1.034482f, 1.041446f, 1.048424f, 1.055417f, 1.062424f,
00173 1.069446f, 1.076482f, 1.083534f, 1.090600f, 1.097681f,
00174 1.104778f, 1.111889f, 1.119016f, 1.126159f, 1.133316f,
00175 1.140490f, 1.147679f, 1.154884f, 1.162105f, 1.169342f,
00176 1.176595f, 1.183864f, 1.191149f, 1.198451f, 1.205770f,
00177 1.213105f, 1.220457f, 1.227826f, 1.235211f, 1.242614f,
00178 1.250034f, 1.257471f, 1.264926f, 1.272398f, 1.279888f,
00179 1.287395f, 1.294921f, 1.302464f, 1.310026f, 1.317605f,
00180 1.325203f, 1.332819f, 1.340454f, 1.348108f, 1.355780f,
00181
00182 1.363472f, 1.371182f, 1.378912f, 1.386660f, 1.394429f,
00183 1.402216f, 1.410024f, 1.417851f, 1.425698f, 1.433565f,
00184 1.441453f, 1.449360f, 1.457288f, 1.465237f, 1.473206f,
00185 1.481196f, 1.489208f, 1.497240f, 1.505293f, 1.513368f,
00186 1.521465f, 1.529583f, 1.537723f, 1.545885f, 1.554068f,
00187 1.562275f, 1.570503f, 1.578754f, 1.587028f, 1.595325f,
00188 1.603644f, 1.611987f, 1.620353f, 1.628743f, 1.637156f,
00189 1.645593f, 1.654053f, 1.662538f, 1.671047f, 1.679581f,
00190 1.688139f, 1.696721f, 1.705329f, 1.713961f, 1.722619f,
00191 1.731303f, 1.740011f, 1.748746f, 1.757506f, 1.766293f,
00192 1.775106f, 1.783945f, 1.792810f, 1.801703f, 1.810623f,
00193 1.819569f, 1.828543f, 1.837545f, 1.846574f, 1.855631f,
00194 1.864717f, 1.873830f, 1.882972f, 1.892143f, 1.901343f,
00195 1.910572f, 1.919830f, 1.929117f, 1.938434f, 1.947781f,
00196 1.957158f, 1.966566f, 1.976004f, 1.985473f, 1.994972f,
00197 2.004503f, 2.014065f, 2.023659f, 2.033285f, 2.042943f,
00198 2.052633f, 2.062355f, 2.072110f, 2.081899f, 2.091720f,
00199 2.101575f, 2.111464f, 2.121386f, 2.131343f, 2.141334f,
00200 2.151360f, 2.161421f, 2.171517f, 2.181648f, 2.191815f,
00201 2.202018f, 2.212257f, 2.222533f, 2.232845f, 2.243195f,
00202
00203 2.253582f, 2.264006f, 2.274468f, 2.284968f, 2.295507f,
00204 2.306084f, 2.316701f, 2.327356f, 2.338051f, 2.348786f,
00205 2.359562f, 2.370377f, 2.381234f, 2.392131f, 2.403070f,
00206 2.414051f, 2.425073f, 2.436138f, 2.447246f, 2.458397f,
00207 2.469591f, 2.480828f, 2.492110f, 2.503436f, 2.514807f,
00208 2.526222f, 2.537684f, 2.549190f, 2.560743f, 2.572343f,
00209 2.583989f, 2.595682f, 2.607423f, 2.619212f, 2.631050f,
00210 2.642936f, 2.654871f, 2.666855f, 2.678890f, 2.690975f,
00211 2.703110f, 2.715297f, 2.727535f, 2.739825f, 2.752168f,
00212 2.764563f, 2.777012f, 2.789514f, 2.802070f, 2.814681f,
00213 2.827347f, 2.840069f, 2.852846f, 2.865680f, 2.878570f,
00214 2.891518f, 2.904524f, 2.917588f, 2.930712f, 2.943894f,
00215 2.957136f, 2.970439f, 2.983802f, 2.997227f, 3.010714f,
00216 3.024263f, 3.037875f, 3.051551f, 3.065290f, 3.079095f,
00217 3.092965f, 3.106900f, 3.120902f, 3.134971f, 3.149107f,
00218 3.163312f, 3.177585f, 3.191928f, 3.206340f, 3.220824f,
00219 3.235378f, 3.250005f, 3.264704f, 3.279477f, 3.294323f,
00220 3.309244f, 3.324240f, 3.339312f, 3.354461f, 3.369687f,
00221 3.384992f, 3.400375f, 3.415838f, 3.431381f, 3.447005f,
00222 3.462711f, 3.478500f, 3.494372f, 3.510328f, 3.526370f,
00223
00224 3.542497f, 3.558711f, 3.575012f, 3.591402f, 3.607881f,
00225 3.624450f, 3.641111f, 3.657863f, 3.674708f, 3.691646f,
00226 3.708680f, 3.725809f, 3.743034f, 3.760357f, 3.777779f,
00227 3.795300f, 3.812921f, 3.830645f, 3.848470f, 3.866400f,
00228 3.884434f, 3.902574f, 3.920821f, 3.939176f, 3.957640f,
00229 3.976215f, 3.994901f, 4.013699f, 4.032612f, 4.051639f,
00230 4.070783f, 4.090045f, 4.109425f, 4.128925f, 4.148547f,
00231 4.168292f, 4.188160f, 4.208154f, 4.228275f, 4.248524f,
00232 4.268903f, 4.289413f, 4.310056f, 4.330832f, 4.351745f,
00233 4.372794f, 4.393982f, 4.415310f, 4.436781f, 4.458395f,
00234 4.480154f, 4.502060f, 4.524114f, 4.546319f, 4.568676f,
00235 4.591187f, 4.613854f, 4.636678f, 4.659662f, 4.682807f,
00236 4.706116f, 4.729590f, 4.753231f, 4.777041f, 4.801024f,
00237 4.825179f, 4.849511f, 4.874020f, 4.898710f, 4.923582f,
00238 4.948639f, 4.973883f, 4.999316f, 5.024942f, 5.050761f,
00239 5.076778f, 5.102993f, 5.129411f, 5.156034f, 5.182864f,
00240 5.209903f, 5.237156f, 5.264625f, 5.292312f, 5.320220f,
00241 5.348354f, 5.376714f, 5.405306f, 5.434131f, 5.463193f,
00242 5.492496f, 5.522042f, 5.551836f, 5.581880f, 5.612178f,
00243 5.642734f, 5.673552f, 5.704634f, 5.735986f, 5.767610f,
00244
00245 5.799512f, 5.831694f, 5.864161f, 5.896918f, 5.929968f,
00246 5.963316f, 5.996967f, 6.030925f, 6.065194f, 6.099780f,
00247 6.134687f, 6.169921f, 6.205486f, 6.241387f, 6.277630f,
00248 6.314220f, 6.351163f, 6.388465f, 6.426130f, 6.464166f,
00249 6.502578f, 6.541371f, 6.580553f, 6.620130f, 6.660109f,
00250 6.700495f, 6.741297f, 6.782520f, 6.824173f, 6.866262f,
00251 6.908795f, 6.951780f, 6.995225f, 7.039137f, 7.083525f,
00252 7.128398f, 7.173764f, 7.219632f, 7.266011f, 7.312910f,
00253 7.360339f, 7.408308f, 7.456827f, 7.505905f, 7.555554f,
00254 7.605785f, 7.656608f, 7.708035f, 7.760077f, 7.812747f,
00255 7.866057f, 7.920019f, 7.974647f, 8.029953f, 8.085952f,
00256 8.142657f, 8.200083f, 8.258245f, 8.317158f, 8.376837f,
00257 8.437300f, 8.498562f, 8.560641f, 8.623554f, 8.687319f,
00258 8.751955f, 8.817481f, 8.883916f, 8.951282f, 9.019600f,
00259 9.088889f, 9.159174f, 9.230477f, 9.302822f, 9.376233f,
00260 9.450735f, 9.526355f, 9.603118f, 9.681054f, 9.760191f,
00261 9.840558f, 9.922186f, 10.005107f, 10.089353f, 10.174959f,
00262 10.261958f, 10.350389f, 10.440287f, 10.531693f, 10.624646f,
00263 10.719188f, 10.815362f, 10.913214f, 11.012789f, 11.114137f,
00264 11.217307f, 11.322352f, 11.429325f, 11.538283f, 11.649285f,
00265
00266 11.762390f, 11.877664f, 11.995170f, 12.114979f, 12.237161f,
00267 12.361791f, 12.488946f, 12.618708f, 12.751161f, 12.886394f,
00268 13.024498f, 13.165570f, 13.309711f, 13.457026f, 13.607625f,
00269 13.761625f, 13.919145f, 14.080314f, 14.245263f, 14.414134f,
00270 14.587072f, 14.764233f, 14.945778f, 15.131877f, 15.322712f,
00271 15.518470f, 15.719353f, 15.925570f, 16.137345f, 16.354912f,
00272 16.578520f, 16.808433f, 17.044929f, 17.288305f, 17.538873f,
00273 17.796967f, 18.062943f, 18.337176f, 18.620068f, 18.912049f,
00274 19.213574f, 19.525133f, 19.847249f, 20.180480f, 20.525429f,
00275 20.882738f, 21.253102f, 21.637266f, 22.036036f, 22.450278f,
00276 22.880933f, 23.329017f, 23.795634f, 24.281981f, 24.789364f,
00277 25.319207f, 25.873062f, 26.452634f, 27.059789f, 27.696581f,
00278 28.365274f, 29.068370f, 29.808638f, 30.589157f, 31.413354f,
00279 32.285060f, 33.208568f, 34.188705f, 35.230920f, 36.341388f,
00280 37.527131f, 38.796172f, 40.157721f, 41.622399f, 43.202525f,
00281 44.912465f, 46.769077f, 48.792279f, 51.005773f, 53.437996f,
00282 56.123356f, 59.103894f,
00283
00284 };
00285
00286 double RandLandau::transform (double r) {
00287
00288 double u = r * TABLE_MULTIPLIER;
00289 int index = int(u);
00290 double du = u - index;
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309 if ( index >= 70 && index <= 800 ) {
00310
00311 double f0 = inverseLandau [index];
00312 double f1 = inverseLandau [index+1];
00313 return f0 + du * (f1 - f0);
00314
00315 } else if ( index >= 7 && index <= 980 ) {
00316
00317 double f_1 = inverseLandau [index-1];
00318 double f0 = inverseLandau [index];
00319 double f1 = inverseLandau [index+1];
00320 double f2 = inverseLandau [index+2];
00321
00322 return f0 + du * (f1 - f0 - .25*(1-du)* (f2 -f1 - f0 + f_1) );
00323
00324 } else if ( index < 7 ) {
00325
00326 const double n0 = 0.99858950;
00327 const double n1 = 34.5213058; const double d1 = 34.1760202;
00328 const double n2 = 17.0854528; const double d2 = 4.01244582;
00329
00330 double logr = std::log(r);
00331 double x = 1/logr;
00332 double x2 = x*x;
00333
00334 double pade = (n0 + n1*x + n2*x2) / (1.0 + d1*x + d2*x2);
00335
00336 return ( - std::log ( -.91893853 - logr ) -1 ) * pade;
00337
00338 } else if ( index <= 999 ) {
00339
00340 const double n0 = 1.00060006;
00341 const double n1 = 263.991156; const double d1 = 257.368075;
00342 const double n2 = 4373.20068; const double d2 = 3414.48018;
00343
00344 double x = 1-r;
00345 double x2 = x*x;
00346
00347 return (n0 + n1*x + n2*x2) / (x * (1.0 + d1*x + d2*x2));
00348
00349 } else {
00350
00351 const double n0 = 1.00001538;
00352 const double n1 = 6075.14119; const double d1 = 6065.11919;
00353 const double n2 = 734266.409; const double d2 = 694021.044;
00354
00355 double x = 1-r;
00356 double x2 = x*x;
00357
00358 return (n0 + n1*x + n2*x2) / (x * (1.0 + d1*x + d2*x2));
00359
00360 }
00361
00362 }
00363
00364 std::ostream & RandLandau::put ( std::ostream & os ) const {
00365 int pr=os.precision(20);
00366 os << " " << name() << "\n";
00367 os.precision(pr);
00368 return os;
00369 }
00370
00371 std::istream & RandLandau::get ( std::istream & is ) {
00372 std::string inName;
00373 is >> inName;
00374 if (inName != name()) {
00375 is.clear(std::ios::badbit | is.rdstate());
00376 std::cerr << "Mismatch when expecting to read state of a "
00377 << name() << " distribution\n"
00378 << "Name found was " << inName
00379 << "\nistream is left in the badbit state\n";
00380 return is;
00381 }
00382 return is;
00383 }
00384
00385 }