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