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