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