GooFit  v2.1.3
VoigtianPdf.cu
Go to the documentation of this file.
1 #include <goofit/Faddeeva.h>
2 #include <goofit/PDFs/basic/VoigtianPdf.h>
3 #include <goofit/detail/Complex.h>
4 #include <limits>
5 
6 namespace GooFit {
7 
8 #define M_2PI 6.28318530717958
9 //#define ROOT2 1.41421356
10 
11 // tables for Pade approximation
12 __constant__ fptype C[7]
13  = {65536.0, -2885792.0, 69973904.0, -791494704.0, 8962513560.0, -32794651890.0, 175685635125.0};
14 __constant__ fptype D[7]
15  = {192192.0, 8648640.0, 183783600.0, 2329725600.0, 18332414100.0, 84329104860.0, 175685635125.0};
16 
17 //#define UNROLL_LOOP 1
18 
19 #ifndef UNROLL_LOOP
20 __constant__ fptype n1[12] = {0.25, 1.0, 2.25, 4.0, 6.25, 9.0, 12.25, 16.0, 20.25, 25.0, 30.25, 36.0};
21 __constant__ fptype e1[12] = {0.7788007830714049,
22  0.3678794411714423,
23  1.053992245618643e-1,
24  1.831563888873418e-2,
25  1.930454136227709e-3,
26  1.234098040866795e-4,
27  4.785117392129009e-6,
28  1.125351747192591e-7,
29  1.605228055185612e-9,
30  1.388794386496402e-11,
31  7.287724095819692e-14,
32  2.319522830243569e-16};
33 
34 // table 2: coefficients for h = 0.53
35 __constant__ fptype n2[12]
36  = {0.2809, 1.1236, 2.5281, 4.4944, 7.0225, 10.1124, 13.7641, 17.9776, 22.7529, 28.09, 33.9889, 40.4496};
37 __constant__ fptype e2[12] = {0.7551038420890235,
38  0.3251072991205958,
39  7.981051630007964e-2,
40  1.117138143353082e-2,
41  0.891593719995219e-3,
42  4.057331392320188e-5,
43  1.052755021528803e-6,
44  1.557498087816203e-8,
45  1.313835773243312e-10,
46  6.319285885175346e-13,
47  1.733038792213266e-15,
48  2.709954036083074e-18};
49 
50 __device__ fpcomplex device_Faddeeva_2(const fpcomplex &z) {
51  fptype *n, *e, t, u, r, s, d, f, g, h;
52  fpcomplex c, d2, v;
53  int i;
54 
55  s = thrust::norm(z); // NB: norm2 is correct although CPU version calls the function 'norm'.
56 
57  if(s < 1e-7) {
58  // use Pade approximation
59  fpcomplex zz = z * z;
60  v = exp(zz); // Note lower-case! This is our own already-templated exp function for thrust::complex, no need for
61  // float/double define.
62  c = C[0];
63  d2 = D[0];
64 
65  for(i = 1; i <= 6; i++) {
66  c = c * zz + C[i];
67  d2 = d2 * zz + D[i];
68  }
69 
70  return fptype(1.0) / v + fpcomplex(0.0, M_2_SQRTPI) * c / d2 * z * v;
71  }
72 
73  // use trapezoid rule
74  // select default table 1
75  n = n1;
76  e = e1;
77  r = M_1_PI * 0.5;
78 
79  // if z is too close to a pole select table 2
80  if(fabs(z.imag()) < 0.01 && fabs(z.real()) < 6.01) {
81  // h = modf(2*fabs(z.real),&g);
82  // Equivalent to above. Do this way because nvcc only knows about double version of modf.
83  h = fabs(z.real()) * 2;
84  g = floor(h);
85  h -= g;
86 
87  if(h < 0.02 || h > 0.98) {
88  n = n2;
89  e = e2;
90  r = M_1_PI * 0.53;
91  }
92  }
93 
94  d = (z.imag() - z.real()) * (z.imag() + z.real());
95  f = 4 * z.real() * z.real() * z.imag() * z.imag();
96 
97  g = h = 0.0;
98 
99  for(i = 0; i < 12; i++) {
100  t = d + n[i];
101  u = e[i] / (t * t + f);
102  g += (s + n[i]) * u;
103  h += (s - n[i]) * u;
104  }
105 
106  u = 1 / s;
107 
108  c = r * fpcomplex(z.imag() * (u + 2.0 * g), z.real() * (u + 2.0 * h));
109 
110  if(z.imag() < M_2PI) {
111  s = 2.0 / r;
112  t = s * z.real();
113  u = s * z.imag();
114  s = sin(t);
115  h = cos(t);
116  f = exp(-u) - h;
117  g = 2.0 * exp(d - u) / (s * s + f * f);
118  u = 2.0 * z.real() * z.imag();
119  h = cos(u);
120  t = sin(u);
121  c += g * fpcomplex((h * f - t * s), -(h * s + t * f));
122  }
123 
124  return c;
125 }
126 
127 #else
128 __device__ fpcomplex device_Faddeeva_2(const fpcomplex &z) {
129  fptype u, s, d, f, g, h;
130  fpcomplex c, d2, v;
131 
132  s = norm2(z); // NB: norm2 is correct although CPU version calls the function 'norm'.
133 
134  if(s < 1e-7) {
135  // use Pade approximation
136  fpcomplex zz = z * z;
137  v = exp(zz); // Note lower-case! This is our own already-templated exp function for thrust::complex, no need for
138  // float/double define.
139  c = C[0];
140  d2 = D[0];
141 
142  for(int i = 1; i < 7; ++i) {
143  c = c * zz + C[i];
144  d2 = d2 * zz + D[i];
145  }
146 
147  return fptype(1.0) / v + fpcomplex(0.0, M_2_SQRTPI) * c / d2 * z * v;
148  }
149 
150  // use trapezoid rule
151  fptype r = M_1_PI * 0.5;
152  bool useDefault = true;
153 
154  // if z is too close to a pole select table 2
155  if(fabs(z.imag) < 0.01 && fabs(z.real) < 6.01) {
156  // h = modf(2*fabs(z.real),&g);
157  // Equivalent to above. Do this way because nvcc only knows about double version of modf.
158  h = fabs(z.real) * 2;
159  g = floor(h);
160  h -= g;
161 
162  if(h < 0.02 || h > 0.98) {
163  useDefault = false;
164  r = M_1_PI * 0.53;
165  }
166  }
167 
168  d = (z.imag - z.real) * (z.imag + z.real);
169  f = 4 * z.real * z.real * z.imag * z.imag;
170 
171  g = h = 0.0;
172  fptype currentN = (useDefault ? 0.25 : 0.2809);
173  fptype currentE = (useDefault ? 0.7788007830714049 : 0.7551038420890235);
174  fptype t = d + currentN;
175  u = currentE / (t * t + f);
176  g += (s + currentN) * u;
177  h += (s - currentN) * u;
178 
179  currentN = (useDefault ? 1.0 : 1.1236);
180  currentE = (useDefault ? 0.3678794411714423 : 0.3251072991205958);
181  t = d + currentN;
182  u = currentE / (t * t + f);
183  g += (s + currentN) * u;
184  h += (s - currentN) * u;
185 
186  currentN = (useDefault ? 2.25 : 2.5281);
187  currentE = (useDefault ? 1.053992245618643e-1 : 7.981051630007964e-2);
188  t = d + currentN;
189  u = currentE / (t * t + f);
190  g += (s + currentN) * u;
191  h += (s - currentN) * u;
192 
193  currentN = (useDefault ? 4.0 : 4.4944);
194  currentE = (useDefault ? 1.930454136227709e-3 : 0.891593719995219e-3);
195  t = d + currentN;
196  u = currentE / (t * t + f);
197  g += (s + currentN) * u;
198  h += (s - currentN) * u;
199 
200  currentN = (useDefault ? 6.25 : 7.0225);
201  currentE = (useDefault ? 4.785117392129009e-6 : 1.052755021528803e-6);
202  t = d + currentN;
203  u = currentE / (t * t + f);
204  g += (s + currentN) * u;
205  h += (s - currentN) * u;
206 
207  currentN = (useDefault ? 9.0 : 10.1124);
208  currentE = (useDefault ? 1.605228055185612e-9 : 1.313835773243312e-10);
209  t = d + currentN;
210  u = currentE / (t * t + f);
211  g += (s + currentN) * u;
212  h += (s - currentN) * u;
213 
214  currentN = (useDefault ? 12.25 : 13.7641);
215  currentE = (useDefault ? 7.287724095819692e-14 : 1.733038792213266e-15);
216  t = d + currentN;
217  u = currentE / (t * t + f);
218  g += (s + currentN) * u;
219  h += (s - currentN) * u;
220 
221  currentN = (useDefault ? 16.0 : 17.9776);
222  currentE = (useDefault ? 1.831563888873418e-2 : 1.117138143353082e-2);
223  t = d + currentN;
224  u = currentE / (t * t + f);
225  g += (s + currentN) * u;
226  h += (s - currentN) * u;
227 
228  currentN = (useDefault ? 20.25 : 22.7529);
229  currentE = (useDefault ? 1.234098040866795e-4 : 4.057331392320188e-5);
230  t = d + currentN;
231  u = currentE / (t * t + f);
232  g += (s + currentN) * u;
233  h += (s - currentN) * u;
234 
235  currentN = (useDefault ? 25.0 : 28.09);
236  currentE = (useDefault ? 1.125351747192591e-7 : 1.557498087816203e-8);
237  t = d + currentN;
238  u = currentE / (t * t + f);
239  g += (s + currentN) * u;
240  h += (s - currentN) * u;
241 
242  currentN = (useDefault ? 30.25 : 33.9889);
243  currentE = (useDefault ? 1.388794386496402e-11 : 6.319285885175346e-13);
244  t = d + currentN;
245  u = currentE / (t * t + f);
246  g += (s + currentN) * u;
247  h += (s - currentN) * u;
248 
249  currentN = (useDefault ? 36.0 : 40.4496);
250  currentE = (useDefault ? 2.319522830243569e-16 : 2.709954036083074e-18);
251  t = d + currentN;
252  u = currentE / (t * t + f);
253  g += (s + currentN) * u;
254  h += (s - currentN) * u;
255 
256  u = 1 / s;
257  c = r * fpcomplex(z.imag * (u + 2.0 * g), z.real * (u + 2.0 * h));
258 
259  if(z.imag < M_2PI) {
260  s = 2.0 / r;
261  t = s * z.real;
262  u = s * z.imag;
263  s = sin(t);
264  h = cos(t);
265  f = exp(-u) - h;
266  g = 2.0 * exp(d - u) / (s * s + f * f);
267  u = 2.0 * z.real * z.imag;
268  h = cos(u);
269  t = sin(u);
270  c += g * fpcomplex((h * f - t * s), -(h * s + t * f));
271  }
272 
273  return c;
274 }
275 #endif
276 
277 __device__ fptype device_Voigtian(fptype *evt, fptype *p, unsigned int *indices) {
278  fptype x = evt[0];
279  fptype m = p[indices[1]];
280  fptype w = p[indices[2]];
281  fptype s = p[indices[3]];
282 
283  // return constant for zero width and sigma
284  if((0 == s) && (0 == w))
285  return 1;
286 
287  // assert(s > 0); // No device-side assert?!
288  // assert(w > 0);
289 
290  fptype arg = x - m;
291 
292  // Breit-Wigner for zero sigma
293  if(0 == s)
294  return (1 / (arg * arg + 0.25 * w * w));
295 
296  fptype coef = -0.5 / (s * s);
297 
298  // Gauss for zero width
299  if(0 == w)
300  return exp(coef * arg * arg);
301 
302  // actual Voigtian for non-trivial width and sigma
303  // fptype c = 1./(ROOT2*s);
304  fptype c = 0.707106781187; // 1/root(2)
305  c /= s;
306  fptype a = 0.5 * c * w;
307  fptype u = c * arg;
308  fpcomplex z(u, a);
309  fpcomplex v = device_Faddeeva_2(z);
310 
311 #define rsqrtPi 0.5641895835477563
312  return c * rsqrtPi * v.real();
313 }
314 
315 __device__ device_function_ptr ptr_to_Voigtian = device_Voigtian;
316 
317 __host__ VoigtianPdf::VoigtianPdf(std::string n, Observable _x, Variable m, Variable s, Variable w)
318  : GooPdf(n, _x) {
319  std::vector<unsigned int> pindices;
320  pindices.push_back(registerParameter(m));
321  pindices.push_back(registerParameter(s));
322  pindices.push_back(registerParameter(w));
323  GET_FUNCTION_ADDR(ptr_to_Voigtian);
324  initialize(pindices);
325 }
326 
327 } // namespace GooFit
Generated on Sat Apr 21 2018 20:13:18 for GooFit by   doxygen 1.8.13