1
2
3
4
5
6
7
8
9
10
11
12 package org.lcsim.math.distribution;
13
14 import static java.lang.Math.sqrt;
15 import static java.lang.Math.exp;
16 import static java.lang.Math.log;
17
18
19
20
21 public class LandauDistribution
22 {
23 private double _peak;
24 private double _width;
25
26
27 public LandauDistribution(double peak, double width)
28 {
29 _peak = peak;
30 _width = width;
31 }
32
33 public double peak()
34 {
35 return _peak;
36 }
37
38 public double width()
39 {
40 return _width;
41 }
42
43 public double value(double x)
44 {
45 double xs = _peak + 0.222782*_width;
46 return denlan((x-xs)/_width)/_width;
47 }
48
49 private double denlan(double x)
50 {
51
52
53
54 double ret_val, r__1;
55
56
57 double u, v;
58 v = x;
59 if (v < -5.5)
60 {
61 u = exp(v + 1.);
62 ret_val = exp(-1 / u) / sqrt(u) * .3989422803 * ((a1[0] + (a1[
63 1] + a1[2] * u) * u) * u + 1);
64 }
65 else if (v < -1.)
66 {
67 u = exp(-v - 1);
68 ret_val = exp(-u) * sqrt(u) * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[
69 4] * v) * v) * v) * v) / (q1[0] + (q1[1] + (q1[2] + (q1[3] +
70 q1[4] * v) * v) * v) * v);
71 }
72 else if (v < 1.)
73 {
74 ret_val = (p2[0] + (p2[1] + (p2[2] + (p2[3] + p2[4] * v) * v) * v) *
75 v) / (q2[0] + (q2[1] + (q2[2] + (q2[3] + q2[4] * v) * v) * v)
76 * v);
77 }
78 else if (v < 5.)
79 {
80 ret_val = (p3[0] + (p3[1] + (p3[2] + (p3[3] + p3[4] * v) * v) * v) *
81 v) / (q3[0] + (q3[1] + (q3[2] + (q3[3] + q3[4] * v) * v) * v)
82 * v);
83 }
84 else if (v < 12.)
85 {
86 u = 1 / v;
87
88 r__1 = u;
89 ret_val = r__1 * r__1 * (p4[0] + (p4[1] + (p4[2] + (p4[3] + p4[4] * u)
90 * u) * u) * u) / (q4[0] + (q4[1] + (q4[2] + (q4[3] + q4[4] *
91 u) * u) * u) * u);
92 }
93 else if (v < 50.)
94 {
95 u = 1 / v;
96
97 r__1 = u;
98 ret_val = r__1 * r__1 * (p5[0] + (p5[1] + (p5[2] + (p5[3] + p5[4] * u)
99 * u) * u) * u) / (q5[0] + (q5[1] + (q5[2] + (q5[3] + q5[4] *
100 u) * u) * u) * u);
101 }
102 else if (v < 300.)
103 {
104 u = 1 / v;
105
106 r__1 = u;
107 ret_val = r__1 * r__1 * (p6[0] + (p6[1] + (p6[2] + (p6[3] + p6[4] * u)
108 * u) * u) * u) / (q6[0] + (q6[1] + (q6[2] + (q6[3] + q6[4] *
109 u) * u) * u) * u);
110 }
111 else
112 {
113 u = 1 / (v - v * log(v) / (v + 1));
114
115 r__1 = u;
116 ret_val = r__1 * r__1 * ((a2[0] + a2[1] * u) * u + 1);
117 }
118 return ret_val;
119
120 }
121
122 static final double[] p1 = { .4259894875,-.124976255,.039842437,-.006298287635,.001511162253 };
123 static final double[] q5 = { 1.,156.9424537,3745.310488,9834.698876,66924.28357 };
124 static final double[] p6 = { 1.000827619,664.9143136,62972.92665,475554.6998,-5743609.109 };
125 static final double[] q6 = { 1.,651.4101098,56974.73333,165917.4725,-2815759.939 };
126 static final double[] a1 = { .04166666667,-.01996527778,.02709538966 };
127 static final double[] a2 = { -1.84556867,-4.284640743 };
128 static final double[] q1 = { 1.,-.3388260629,.09594393323,-.01608042283,.003778942063 };
129 static final double[] p2 = { .1788541609,.1173957403,.01488850518,-.001394989411,1.283617211e-4 };
130 static final double[] q2 = { 1.,.7428795082,.3153932961,.06694219548,.008790609714 };
131 static final double[] p3 = { .1788544503,.09359161662,.006325387654,6.611667319e-5,-2.031049101e-6 };
132 static final double[] q3 = { 1.,.6097809921,.2560616665,.04746722384,.006957301675 };
133 static final double[] p4 = { .9874054407,118.6723273,849.279436,-743.7792444,427.0262186 };
134 static final double[] q4 = { 1.,106.8615961,337.6496214,2016.712389,1597.063511 };
135 static final double[] p5 = { 1.003675074,167.5702434,4789.711289,21217.86767,-22324.9491 };
136
137 }