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