GooFit  v2.1.3
LineshapesPdf.cu
Go to the documentation of this file.
1 /*
2 04/05/2016 Christoph Hasse
3 DISCLAIMER:
4 
5 This code is not sufficently tested yet and still under heavy development!
6 
7 This file includes some lineshapes and spinfactors.
8 Also right now it is the home to some helper functions needed and an implementation of a simple 4-vec class that works
9 on the GPU
10 */
11 
12 #include <goofit/PDFs/physics/LineshapesPdf.h>
13 #include <goofit/PDFs/physics/SpinFactors.h>
14 #include <goofit/Version.h>
15 
16 #if GOOFIT_KMATRIX && THRUST_DEVICE_SYSTEM == THRUST_DEVICE_SYSTEM_CUDA
17 #include <goofit/detail/compute_inverse5.h>
18 #endif
19 
20 #include <utility>
21 
22 #include <Eigen/Core>
23 #include <Eigen/LU>
24 
25 #include <goofit/detail/Macros.h>
26 
27 #define NPOLES 5
28 #define NCHANNELS 5
29 
30 #define mPiPlus 0.139570
31 #define mKPlus 0.493677
32 #define mEta 0.547862
33 #define mEtap 0.96778
34 
35 namespace GooFit {
36 
37 // Lineshape base
38 
39 // Form factors as in pdg http://pdg.lbl.gov/2012/reviews/rpp2012-rev-dalitz-analysis-formalism.pdf
40 __device__ fptype BL_PRIME(fptype z2, fptype z02, int L) {
41  if(0 == L)
42  return 1.0;
43  else if(1 == L)
44  return (1 + z02) / (1 + z2);
45  else if(2 == L)
46  return (z02 * z02 + 3 * z02 + 9) / (z2 * z2 + 3 * z2 + 9);
47  else {
48  printf("ERROR! Oribtal > 2 not supported!\n");
49  return 0;
50  }
51 
52  // Spin 3 and up not accounted for.
53 }
54 
55 __device__ fptype BL(fptype z2, int L) {
56  if(0 == L)
57  return 1.0;
58  else if(1 == L)
59  return 2 * z2 / (1 + z2);
60  else if(2 == L)
61  return (13 * z2 * z2) / (z2 * z2 + 3 * z2 + 9);
62  else {
63  printf("ERROR! Oribtal > 2 not supported!\n");
64  return 0;
65  }
66 
67  // Spin 3 and up not accounted for.
68 }
69 
70 __device__ fptype BL2(fptype z2, int L) {
71  if(0 == L)
72  return 1.0;
73  else if(1 == L)
74  return 1.0 / (1 + z2);
75  else if(2 == L)
76  return 1.0 / (z2 * z2 + 3 * z2 + 9);
77  else {
78  printf("ERROR! Oribtal > 2 not supported!\n");
79  return 0;
80  }
81 
82  // Spin 3 and up not accounted for.
83 }
84 
85 __device__ fpcomplex LS_ONE(fptype Mpair, fptype m1, fptype m2, unsigned int *indices) { return {1., 0.}; }
86 
87 // This function is modeled after BW_BW::getVal() in BW_BW.cpp from the MINT package written by Jonas Rademacker.
88 __device__ fpcomplex BW(fptype Mpair, fptype m1, fptype m2, unsigned int *indices) {
89  fptype meson_radius = functorConstants[indices[7]];
90  fptype resmass = cudaArray[indices[2]];
91  fptype reswidth = cudaArray[indices[3]];
92  unsigned int orbital = indices[4];
93  unsigned int FF = indices[6];
94 
95  const unsigned int to2Lplus1 = 2 * orbital + 1;
96 
97  fptype mass = resmass;
98  fptype width = reswidth;
99  fptype mumsRecoMass2 = Mpair * Mpair;
100 
101  fptype mpsq = (m1 + m2) * (m1 + m2);
102  fptype mmsq = (m1 - m2) * (m1 - m2);
103  fptype num = (mumsRecoMass2 - mpsq) * (mumsRecoMass2 - mmsq);
104  fptype num2 = (mass * mass - mpsq) * (mass * mass - mmsq);
105  fptype pABSq = num / (4 * mumsRecoMass2);
106  fptype prSqForGofM = num2 / (4 * mass * mass);
107  fptype prSq2 = prSqForGofM < 0 ? 0 : prSqForGofM;
108  prSqForGofM = fabs(prSqForGofM);
109 
110  fptype pratio = sqrt(pABSq / prSqForGofM);
111 
112  fptype pratio_to_2Jplus1 = 1;
113 
114  for(int i = 0; i < to2Lplus1; i++) {
115  pratio_to_2Jplus1 *= pratio;
116  }
117 
118  fptype mratio = mass / Mpair;
119  fptype r2 = meson_radius * meson_radius;
120  fptype thisFR = BL_PRIME(pABSq * r2, prSqForGofM * r2, orbital);
121  fptype frFactor = 1;
122 
123  if(0 != orbital && 0 != FF) {
124  frFactor = (FF == 1 ? BL(pABSq * r2, orbital) : BL_PRIME(pABSq * r2, prSq2 * r2, orbital));
125  frFactor = (FF == 3 ? BL2(pABSq * r2, orbital) : frFactor);
126  }
127 
128  fptype GofM = width * pratio_to_2Jplus1 * mratio * thisFR;
129 
130  fptype gamma = mass * sqrt((mass * mass + width * width));
131  fptype k = (2.0 * sqrt(2.0) / M_PI) * mass * width * gamma
132  / sqrt(mass * mass + gamma); // Note added additional factor of 2*sqrt(2)/PI here so results are
133  // comparable to MINT3. MINT2 doesn't have include this.
134 
135  fpcomplex BW(mass * mass - mumsRecoMass2, mass * GofM);
136  fptype den = (mass * mass - mumsRecoMass2) * (mass * mass - mumsRecoMass2) + mass * GofM * mass * GofM;
137 
138  fpcomplex ret = (sqrt(k * frFactor)) / den * BW;
139  // printf("m1, m2, Mpair, to2Lplus1, GofM, thisFR, pratio, mratio, pABSq , prSqForGofM, FF, ret.real, ret.imag\n");
140  // printf("BW %.7g, %.7g, %.7g, %i, %i, %i, %i\n",meson_radius, resmass, reswidth, orbital, FF, indices[2],
141  // indices[3]);
142  // printf("BW %.7g, %.7g, %.7g, %i, %.7g, %.7g, %.7g, %.7g, %.7g, %.7g, %.7g, %.7g, %.7g\n", m1, m2, Mpair,
143  // to2Lplus1, GofM, thisFR, pratio, mratio, pABSq, prSqForGofM, frFactor, ret.real, ret.imag );
144  return ret;
145 }
146 
147 // This function is modeled after SBW from the MINT package written by Jonas Rademacker.
148 __device__ fpcomplex SBW(fptype Mpair, fptype m1, fptype m2, unsigned int *indices) {
149  fptype resmass = GOOFIT_GET_PARAM(2);
150  fptype reswidth = GOOFIT_GET_PARAM(3);
151  unsigned int orbital = GOOFIT_GET_INT(4);
152  // GOOFIT_GET_INT(5, Mpair, "Mpair");
153  unsigned int FF = GOOFIT_GET_INT(6);
154  fptype meson_radius = GOOFIT_GET_CONST(7);
155 
156  // fptype meson_radius = functorConstants[indices[7]];
157  // fptype resmass = cudaArray[indices[2]];
158  // fptype reswidth = cudaArray[indices[3]];
159  // unsigned int orbital = indices[4];
160  // unsigned int FF = indices[6];
161 
162  fptype mass = resmass;
163  fptype width = reswidth;
164  fptype mumsRecoMass2 = Mpair * Mpair;
165 
166  fptype mpsq = (m1 + m2) * (m1 + m2);
167  fptype mmsq = (m1 - m2) * (m1 - m2);
168  fptype num = (mumsRecoMass2 - mpsq) * (mumsRecoMass2 - mmsq);
169  fptype num2 = (mass * mass - mpsq) * (mass * mass - mmsq);
170  fptype pABSq = num / (4 * mumsRecoMass2);
171  fptype prSq = num2 / (4 * mass * mass);
172  fptype prSq2 = prSq < 0 ? 0 : prSq;
173  prSq = fabs(prSq);
174 
175  fptype r2 = meson_radius * meson_radius;
176  fptype frFactor = 1;
177 
178  if(0 != orbital && 0 != FF) {
179  frFactor = (FF == 1 ? BL(pABSq * r2, orbital) : BL_PRIME(pABSq * r2, prSq2 * r2, orbital));
180  frFactor = (FF == 3 ? BL2(pABSq * r2, orbital) : frFactor);
181  }
182 
183  fptype GofM = width;
184 
185  fptype gamma = sqrt(mass * mass * (mass * mass + width * width));
186  fptype k = mass * width * gamma / sqrt(mass * mass + gamma);
187 
188  fpcomplex BW(mass * mass - mumsRecoMass2, mass * GofM);
189  fptype den = (mass * mass - mumsRecoMass2) * (mass * mass - mumsRecoMass2) + mass * GofM * mass * GofM;
190 
191  fpcomplex ret = (sqrt(k * frFactor)) / den * BW;
192 
193  // printf("m1, m2, Mpair, GofM, pABSq , prSq, FF, ret.real, ret.imag\n");
194  // printf("SBW %.7g, %.7g, %.7g, %.7g, %.7g, %.7g, %.7g, %.7g, %.7g\n", m1, m2, Mpair, GofM, pABSq, prSq, frFactor,
195  // ret.real, ret.imag );
196  return ret;
197 }
198 
199 __device__ fpcomplex bugg_rho2(const fptype &s, const fptype m) {
200  fptype rho_squared = 1. - 4. * m * m / s;
201  fpcomplex returnVal = (rho_squared >= 0) ? fpcomplex(1, 0) : fpcomplex(0, 1);
202  rho_squared = (rho_squared >= 0) ? sqrt(rho_squared) : sqrt(-rho_squared);
203  return rho_squared * returnVal;
204 }
205 
206 __device__ fptype bugg_j1(const fptype &s, const fptype m) {
207  fptype rho_pipi = bugg_rho2(s, m).real();
208  fptype returnVal = 2.;
209  returnVal += (rho_pipi > 0.) ? rho_pipi * log((1. - rho_pipi) / (1. + rho_pipi)) : 0;
210  return returnVal / M_PI;
211 }
212 
213 __device__ fptype bugg_Gamma_4pi(const fptype &s,
214  const fptype mpi,
215  const fptype &g_4pi,
216  const fptype &M,
217  const fptype &lambda_4pi,
218  const fptype &s0_4pi) {
219  fptype returnVal = (s < (16. * mpi * mpi)) ? 0
220  : g_4pi * (1. / (1 + exp(lambda_4pi * (s0_4pi - s))))
221  / (1. / (1 + exp(lambda_4pi * (s0_4pi - M * M))));
222  return returnVal;
223 }
224 
225 // This function is an adaptation from the bugg lineshape implemented in the MINT package written by Jonas Rademacker.
226 // this lineshape is not tested yet!
227 __device__ fpcomplex bugg_MINT(fptype Mpair, fptype m1, fptype m2, unsigned int *indices) {
228  fptype s = Mpair * Mpair;
229 
230  fptype M = 0.953;
231  fptype b1 = 1.302;
232  fptype b2 = 0.340;
233  fptype A = 2.426;
234  fptype g_4pi = 0.011;
235  fptype g_2K = 0.6;
236  fptype g_2eta = 0.2;
237  fptype alpha = 1.3;
238  fptype sA = 0.41;
239  fptype s0_4pi = 7.082 / 2.845;
240  fptype lambda_4pi = 2.845;
241 
242  fptype g1sq = (b1 + b2 * s) * exp(-(s - M * M) / A);
243  fptype z = bugg_j1(s, mPiPlus) - bugg_j1(M * M, mPiPlus);
244 
245  fpcomplex gamma_2pi = fpcomplex(
246  g1sq * (s - sA * mPiPlus * mPiPlus) / (M * M - sA * mPiPlus * mPiPlus) * bugg_rho2(s, mPiPlus).real(), 0);
247  fpcomplex gamma_2K = g_2K * g1sq * s / (M * M)
248  * exp((-1) * alpha * sqrt((s - 4. * mKPlus * mKPlus) * (s - 4. * mKPlus * mKPlus)))
249  * bugg_rho2(s, mKPlus);
250  fpcomplex gamma_2eta = g_2eta * g1sq * s / (M * M)
251  * exp((-1) * alpha * sqrt((s - 4. * mEta * mEta) * (s - 4. * mEta * mEta)))
252  * bugg_rho2(s, mEta);
253  fpcomplex gamma_4pi = fpcomplex(bugg_Gamma_4pi(s, mPiPlus, g_4pi, M, lambda_4pi, s0_4pi), 0);
254 
255  fpcomplex Gamma_tot = gamma_2pi + gamma_2K + gamma_2eta + gamma_4pi;
256 
257  // fpcomplex num = M * gamma_2pi; //only for elastic scattering, not production
258  fpcomplex den
259  = fpcomplex(M * M - s - M * g1sq * (s - sA * mPiPlus * mPiPlus) / (M * M - sA * mPiPlus * mPiPlus) * z, 0)
260  - fpcomplex(0, 1) * M * Gamma_tot;
261  fpcomplex returnVal = 1.0 / den;
262  // printf("Bugg %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g \n",gamma_2pi.real, gamma_2pi.imag, gamma_2K.real,
263  // gamma_2K.imag, gamma_2eta.real, gamma_2eta.imag, gamma_4pi.real, gamma_4pi.imag);
264  // printf("Bugg %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g \n",Mpair, Gamma_tot.real, Gamma_tot.imag, g1sq, z,
265  // den.real, den.imag, returnVal.real, returnVal.imag);
266 
267  // the factor sqrt(1000) gives the correct result in comparison with mint2, I think its because BW/SBW
268  // have a factor of sqrt(k) which these lineshapes dont have. For now this stays because it works. further
269  // investigation needed.
270  return returnVal * sqrt(1000.0);
271 }
272 
273 __device__ fpcomplex bugg_MINT3(fptype Mpair, fptype m1, fptype m2, unsigned int *indices) {
274  fptype s = Mpair * Mpair;
275  fptype M = 0.953;
276  fptype b1 = 1.302;
277  fptype b2 = 0.340;
278  fptype A = 2.426;
279  fptype g_4pi = 0.011;
280  fptype g_2K = 0.6;
281  fptype g_2eta = 0.2;
282  fptype alpha = 1.3;
283  fptype s0_4pi = 7.082 / 2.845;
284  fptype lambda_4pi = 2.845;
285  fptype sA = 0.41 * mPiPlus * mPiPlus;
286 
287  fptype g1sq = M * (b1 + b2 * s) * exp(-(s - M * M) / A);
288  fptype z = bugg_j1(s, mPiPlus) - bugg_j1(M * M, mPiPlus);
289  fptype adlerZero = (s - sA) / (M * M - sA);
290 
291  fptype mk4 = 4. * mKPlus * mKPlus;
292  fptype me4 = 4. * mEta * mEta;
293  fptype tmp1 = s > mk4 ? s - mk4 : mk4 - s;
294  fptype tmp2 = s > me4 ? s - me4 : me4 - s;
295 
296  fpcomplex gamma_2pi = fpcomplex(g1sq * adlerZero * bugg_rho2(s, mPiPlus).real(), 0);
297  fpcomplex gamma_2K = g_2K * g1sq * s / (M * M) * exp((-1) * alpha * tmp1) * bugg_rho2(s, mKPlus);
298  fpcomplex gamma_2eta = g_2eta * g1sq * s / (M * M) * exp((-1) * alpha * tmp2) * bugg_rho2(s, mEta);
299  fpcomplex gamma_4pi = fpcomplex(bugg_Gamma_4pi(s, mPiPlus, g_4pi, M, lambda_4pi, s0_4pi), 0);
300 
301  fpcomplex Gamma_tot = gamma_2pi + gamma_2K + gamma_2eta + gamma_4pi;
302 
303  // fpcomplex num = M * gamma_2pi; //only for elastic scattering, not production
304  fpcomplex den = fpcomplex(M * M - s - adlerZero * g1sq * z, 0) - fpcomplex(0, 1) * Gamma_tot;
305  fpcomplex returnVal = 1.0 / den;
306  // printf("Bugg %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g \n",gamma_2pi.real, gamma_2pi.imag, gamma_2K.real,
307  // gamma_2K.imag, gamma_2eta.real, gamma_2eta.imag, gamma_4pi.real, gamma_4pi.imag);
308  // printf("Bugg %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g \n",Mpair, Gamma_tot.real, Gamma_tot.imag, g1sq, z,
309  // den.real, den.imag, returnVal.real, returnVal.imag);
310 
311  return returnVal;
312 }
313 
314 __device__ fpcomplex lass_MINT(fptype Mpair, fptype m1, fptype m2, unsigned int *indices) {
315  fptype resmass = cudaArray[indices[2]];
316  fptype reswidth = cudaArray[indices[3]];
317  fptype rMass2 = Mpair * Mpair;
318 
319  fptype a = 2.07;
320  fptype r = 3.32;
321 
322  fptype mpsq = (m1 + m2) * (m1 + m2);
323  fptype mmsq = (m1 - m2) * (m1 - m2);
324  fptype num = (rMass2 - mpsq) * (rMass2 - mmsq);
325  fptype num2 = (resmass * resmass - mpsq) * (resmass * resmass - mmsq);
326  fptype pABSq = num / (4 * rMass2);
327  fptype prSq = fabs(num2 / (4 * resmass * resmass));
328 
329  fptype y = 2.0 * a * sqrt(pABSq);
330  fptype x = 2.0 + a * r * pABSq;
331  fptype cotDeltaBg = x / y;
332  fpcomplex phaseshift((cotDeltaBg * cotDeltaBg - 1) / (1 + cotDeltaBg * cotDeltaBg),
333  2 * cotDeltaBg / (1 + cotDeltaBg * cotDeltaBg));
334  // (cotDeltaBg*cotDeltaBg-1)/(1+cotDeltaBg*cotDeltaBg) = cos(2*delta) 2*cotDeltaBg / ( 1 +
335  // cotDeltaBg*cotDeltaBg) = sin(2*delta)
336  fpcomplex den(sqrt(pABSq) * cotDeltaBg, (-1.) * sqrt(pABSq));
337  fptype SF = Mpair * sqrt(prSq) / (resmass * resmass * reswidth);
338  fpcomplex BG = SF / den;
339  fpcomplex returnVal = BG + phaseshift * BW(Mpair, m1, m2, indices);
340  // printf("Lass: %.5g %.5g %.5g %.5g %.5g %.5g\n",BG.real, BG.imag, phaseshift.real, phaseshift.imag,
341  // returnVal.real, returnVal.imag);
342 
343  return returnVal;
344 }
345 
346 // generalized lass lineshape as implemented in MINT3 by Tim Evans. if F=R=1 and phiF=phiR=0 this is equal to normal
347 // lass as implemented in Mint3.
348 // The difference between this and lass mint is not quite clear to me. need to get back to this later.
349 __device__ fpcomplex glass_MINT3(fptype Mpair, fptype m1, fptype m2, unsigned int *indices) {
350  fptype meson_radius = functorConstants[indices[7]];
351  fptype resmass = cudaArray[indices[2]];
352  fptype reswidth = cudaArray[indices[3]];
353  unsigned int orbital = indices[4];
354  fptype rMass2 = Mpair * Mpair;
355 
356  // fptype a = 2.07;
357  // fptype r = 3.32;
358  // fptype phiF = 0.0;
359  // fptype phiR = 0.0;
360  // fptype F = 1.0;
361  fptype a = cudaArray[indices[8]];
362  fptype r = cudaArray[indices[9]];
363  fptype phiF = cudaArray[indices[10]];
364  fptype phiR = cudaArray[indices[11]];
365  fptype F = cudaArray[indices[12]];
366 
367  fptype R = 1.0;
368  // printf("GLass: %.5g %.5g %.5g %.5g %.5g %.5g\n",a, r, phiF, phiR, F, R);
369  // printf("GLass2: %.5g %.5g %.5g %u \n",meson_radius, resmass, reswidth, orbital);
370 
371  fptype mpsq = (m1 + m2) * (m1 + m2);
372  fptype mmsq = (m1 - m2) * (m1 - m2);
373  fptype num = (rMass2 - mpsq) * (rMass2 - mmsq);
374  fptype num2 = (resmass * resmass - mpsq) * (resmass * resmass - mmsq);
375  fptype pABSq = num / (4 * rMass2);
376  fptype prSq = fabs(num2 / (4 * resmass * resmass));
377 
378  fptype pratio = sqrt(pABSq / prSq);
379 
380  fptype pratio_to_2Jplus1 = 1;
381 
382  for(int i = 0; i < 2 * orbital + 1; i++) {
383  pratio_to_2Jplus1 *= pratio;
384  }
385 
386  fptype mratio = resmass / Mpair;
387  fptype r2 = meson_radius * meson_radius;
388  fptype thisFR = BL_PRIME(pABSq * r2, prSq * r2, orbital);
389  fptype GofM = reswidth * pratio_to_2Jplus1 * mratio * thisFR;
390 
391  fptype y = 2.0 * a * sqrt(pABSq);
392  fptype x = 2.0 + a * r * pABSq;
393  fptype scattphase = phiF + atan(y / x);
394  fptype resphase = phiR + atan(resmass * GofM / (resmass * resmass - rMass2));
395  fptype rho = 1.0 / sqrt(pABSq / rMass2);
396  fpcomplex returnVal
397  = (F * sin(scattphase) * fpcomplex(cos(scattphase), sin(scattphase))
398  + R * sin(resphase) * fpcomplex(cos(resphase + 2 * scattphase), sin(resphase + 2 * scattphase)))
399  * rho;
400  // printf("GLass3: %.5g %.5g %.5g %.5g %.5g %.5g\n",rMass2, pABSq, rho, GofM, scattphase, resphase);
401 
402  // printf("GLass4: %.5g %.5g\n",returnVal.real, returnVal.imag);
403  return returnVal;
404 }
405 
406 __device__ fpcomplex aSqrtTerm(const fptype &m0, const fptype &m) {
407  fptype a2 = 1 - (2 * m0 / m) * (2 * m0 / m);
408  fpcomplex returnVal = a2 > 0 ? fpcomplex(sqrt(a2), 0) : fpcomplex(0, sqrt(-a2));
409  return returnVal;
410 }
411 
412 __device__ fpcomplex Flatte_MINT(fptype Mpair, fptype m1, fptype m2, unsigned int *indices) {
413  fptype meson_radius = functorConstants[indices[7]];
414  fptype resmass = cudaArray[indices[2]];
415  unsigned int orbital = indices[4];
416  fptype frFactor = 1;
417  fptype rMass2 = Mpair * Mpair;
418 
419  // As far as I understand, this is only valid for the f980
420  fptype gPi = .165;
421  fptype gK_by_gPi = 4.21;
422  fptype gK = gPi * gK_by_gPi;
423  fptype mPi0 = .1349766;
424  fptype mK0 = .497648;
425 
426  fptype mpsq = (m1 + m2) * (m1 + m2);
427  fptype mmsq = (m1 - m2) * (m1 - m2);
428  fptype num = (rMass2 - mpsq) * (rMass2 - mmsq);
429  // fptype num2 = (resmass*resmass - mpsq)*(resmass*resmass - mmsq);
430  fptype pABSq = num / (4 * rMass2);
431  // fptype prSq = fabs(num2/(4*resmass*resmass));
432 
433  fpcomplex Gpipi = (1. / 3.) * aSqrtTerm(mPi0, Mpair) + (2. / 3.) * aSqrtTerm(mPiPlus, Mpair);
434  fpcomplex GKK = (1. / 2.) * aSqrtTerm(mK0, Mpair) + (1. / 2.) * aSqrtTerm(mKPlus, Mpair);
435  fpcomplex FlatteWidth = gPi * Gpipi + gK * GKK;
436  // printf("%.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g \n",Gpipi.real, Gpipi.imag, GKK.real, GKK.imag, FlatteWidth.real,
437  // FlatteWidth.imag, Mpair, pABSq);
438 
439  frFactor = BL2(pABSq * meson_radius * meson_radius, orbital);
440  fpcomplex BW = sqrt(frFactor) / fpcomplex(resmass * resmass - rMass2, 0) - fpcomplex(0, 1) * resmass * FlatteWidth;
441  return BW;
442 }
443 
444 __device__ __thrust_forceinline__ fptype Q2(fptype Msq, fptype M1sq, fptype M2sq) {
445  return (Msq / 4. - (M1sq + M2sq) / 2. + (M1sq - M2sq) * (M1sq - M2sq) / (4 * Msq));
446 }
447 
448 __device__ __thrust_forceinline__ fptype BlattWeisskopf_Norm(const fptype z2, const fptype z02, unsigned int L) {
449  if(L == 0)
450  return 1;
451  else if(L == 1)
452  return (1 + z02) / (1 + z2);
453  else if(L == 2)
454  return (z02 * z02 + 3 * z02 + 9) / (z2 * z2 + 3 * z2 + 9);
455  else {
456  abort(__FILE__, __LINE__, "Wrong value of L");
457  return 0; // Can't reach
458  }
459 }
460 
461 __device__ fptype getSpline(fptype x, bool continued, unsigned int *indices) {
462  const fptype s_min = GOOFIT_GET_CONST(8);
463  const fptype s_max = GOOFIT_GET_CONST(9);
464  const unsigned int nBins = GOOFIT_GET_INT(10);
465 
466  // 11 is the first spine knot, 11+nBins is the first curvature
467 
468  if(x <= s_min)
469  return continued ? GOOFIT_GET_PARAM(11 + 0) : 0;
470  if(x >= s_max)
471  return continued ? GOOFIT_GET_PARAM(11 + nBins - 1) : 0;
472 
473  fptype spacing = (s_max - s_min) / (nBins - 1.);
474  fptype dx = fmod((x - s_min), spacing);
475 
476  auto bin = static_cast<unsigned int>((x - s_min) / spacing);
477 
478  fptype m_x_0 = GOOFIT_GET_PARAM(11 + bin);
479  fptype m_x_1 = GOOFIT_GET_PARAM(11 + bin + 1);
480  fptype m_xf_0 = GOOFIT_GET_CONST(11 + bin + nBins);
481  fptype m_xf_1 = GOOFIT_GET_CONST(11 + bin + nBins + 1);
482 
483  return m_x_0 + dx * ((m_x_1 - m_x_0) / spacing - (m_xf_1 + 2 * m_xf_0) * spacing / 6) + dx * dx * m_xf_0
484  + dx * dx * dx * (m_xf_1 - m_xf_0) / (6 * spacing);
485 }
486 
487 __device__ fptype kFactor(fptype mass, fptype width) {
488  fptype gamma = mass * sqrt(POW2(mass) + POW2(width));
489  fptype k = 2 * sqrt(2.) * mass * width * gamma / (M_PI * sqrt(POW2(mass) + gamma));
490  return sqrt(k);
491 }
492 
493 __device__ fpcomplex Spline_TDP(fptype Mpair, fptype m1, fptype m2, unsigned int *indices) {
494  const fptype mass = GOOFIT_GET_PARAM(2);
495  const fptype width = GOOFIT_GET_PARAM(3);
496  // const unsigned int L = GOOFIT_GET_INT(4);
497  const fptype radius = GOOFIT_GET_CONST(7);
498 
499  fptype s = POW2(Mpair);
500  fptype s1 = POW2(m1);
501  fptype s2 = POW2(m2);
502 
503  // This is GSpline.EFF in AmpGen
504 
505  fptype q2 = fabs(Q2(s, s1, s2));
506 
507  // Non-EFF
508  // fptype BF = sqrt( BlattWeisskopf_Norm(q2 * POW2(radius), 0, L));
509  fptype BF = exp(-q2 * POW2(radius) / 2);
510 
511  fptype width_shape = width * getSpline(s, true, indices);
512  fptype width_norm = width * getSpline(POW2(mass), false, indices);
513 
514  fptype norm = kFactor(mass, width) * BF;
515  fptype running_width = width * width_shape / width_norm;
516  fpcomplex iBW = fpcomplex(POW2(mass) - s, -mass * running_width);
517  return norm / iBW;
518 }
519 
520 __device__ fpcomplex nonres_DP(fptype Mpair, fptype m1, fptype m2, unsigned int *indices) {
521  fptype meson_radius = functorConstants[indices[7]];
522  unsigned int orbital = indices[4];
523 
524  fptype mumsRecoMass2 = Mpair * Mpair;
525 
526  fptype mpsq = (m1 + m2) * (m1 + m2);
527  fptype mmsq = (m1 - m2) * (m1 - m2);
528  fptype num = (mumsRecoMass2 - mpsq) * (mumsRecoMass2 - mmsq);
529  fptype pABSq = num / (4 * mumsRecoMass2);
530  fptype formfactor = sqrt(BL2(pABSq * meson_radius * meson_radius, orbital));
531  // printf("NonRes q2:%.7g FF:%.7g, s %.7g m1 %.7g m2 %.7g r %.7g L %u \n",pABSq, formfactor, mumsRecoMass2,
532  // m1,m2,meson_radius, orbital );
533  return fpcomplex(1., 0.) * formfactor;
534 }
535 
536 __device__ fptype phsp_twoBody(fptype s, fptype m0, fptype m1) { return sqrt(1. - POW2(m0 + m1) / s); }
537 
538 __device__ fptype phsp_fourPi(fptype s) {
539  if(s > 1)
540  return phsp_twoBody(s, 2 * mPiPlus, 2 * mPiPlus);
541  else
542  return 0.00051 + -0.01933 * s + 0.13851 * s * s + -0.20840 * s * s * s + -0.29744 * s * s * s * s
543  + 0.13655 * s * s * s * s * s + 1.07885 * s * s * s * s * s * s;
544 }
545 
546 #if GOOFIT_KMATRIX
547 
548 __device__ Eigen::Array<fpcomplex, NCHANNELS, NCHANNELS>
549 getPropagator(const Eigen::Array<fptype, NCHANNELS, NCHANNELS> &kMatrix,
550  const Eigen::Matrix<fptype, 5, 1> &phaseSpace,
551  fptype adlerTerm) {
552  Eigen::Array<fpcomplex, NCHANNELS, NCHANNELS> tMatrix;
553 
554  for(unsigned int i = 0; i < NCHANNELS; ++i) {
555  for(unsigned int j = 0; j < NCHANNELS; ++j) {
556  tMatrix(i, j) = (i == j ? 1. : 0.) - fpcomplex(0, adlerTerm) * kMatrix(i, j) * phaseSpace(j);
557  }
558  }
559 
560 #if THRUST_DEVICE_SYSTEM == THRUST_DEVICE_SYSTEM_CUDA
561  // Here we assume that some values are 0
562  return compute_inverse5<-1,
563  -1,
564  0,
565  -1,
566  -1,
567  -1,
568  -1,
569  0,
570  -1,
571  -1,
572  -1,
573  -1,
574  -1,
575  -1,
576  -1,
577  -1,
578  -1,
579  -1,
580  -1,
581  -1,
582  -1,
583  -1,
584  -1,
585  -1,
586  -1>(tMatrix);
587 #else
588  return Eigen::inverse(tMatrix);
589 #endif
590 }
591 
592 __device__ fpcomplex kMatrixFunction(fptype Mpair, fptype m1, fptype m2, unsigned int *indices) {
593  // const fptype mass = GOOFIT_GET_PARAM(2);
594  // const fptype width = GOOFIT_GET_PARAM(3);
595  // const unsigned int L = GOOFIT_GET_INT(4);
596  // const fptype radius = GOOFIT_GET_CONST(7);
597 
598  // const fptype pTerm = GOOFIT_GET_INT();
599 
600  unsigned int pterm = GOOFIT_GET_INT(8);
601  bool is_pole = GOOFIT_GET_INT(9) == 1;
602 
603  fptype sA0 = GOOFIT_GET_PARAM(10);
604  fptype sA = GOOFIT_GET_PARAM(11);
605  fptype s0_prod = GOOFIT_GET_PARAM(12);
606  fptype s0_scatt = GOOFIT_GET_PARAM(13);
607 
608  Eigen::Array<fptype, NCHANNELS, 1> fscat;
609  Eigen::Array<fptype, NPOLES, 1> pmasses;
610  Eigen::Array<fptype, NPOLES, NPOLES> couplings;
611 
612  for(int i = 0; i < NCHANNELS; i++) {
613  fscat(i) = GOOFIT_GET_PARAM(14 + i);
614  }
615 
616  for(int i = 0; i < NPOLES; i++) {
617  for(int j = 0; j < NPOLES; j++)
618  couplings(i, j) = GOOFIT_GET_PARAM(14 + NCHANNELS + i * (NPOLES + 1) + j);
619  pmasses(i) = GOOFIT_GET_PARAM(14 + NCHANNELS + i * (NPOLES + 1) + NPOLES);
620  }
621 
622  fptype s = POW2(Mpair);
623 
624  // constructKMatrix
625 
626  Eigen::Array<fptype, NCHANNELS, NCHANNELS> kMatrix;
627  kMatrix.setZero();
628 
629  // TODO: Make sure the order (k,i,j) is correct
630 
631  for(int i = 0; i < 5; i++) {
632  for(int j = 0; j < 5; j++) {
633  for(int k = 0; k < 5; k++)
634  kMatrix(i, j) += couplings(k, i) * couplings(k, j) / (pmasses(k) - s);
635  if(i == 0 || j == 0) // Scattering term
636  kMatrix(i, j) += fscat(i + j) * (1 - s0_scatt) / (s - s0_scatt);
637  }
638  }
639 
640  fptype adlerTerm = (1. - sA0) * (s - sA * mPiPlus * mPiPlus / 2) / (s - sA0);
641 
642  Eigen::Matrix<fptype, 5, 1> phaseSpace;
643  phaseSpace << phsp_twoBody(s, mPiPlus, mPiPlus), phsp_twoBody(s, mKPlus, mKPlus), phsp_fourPi(s),
644  phsp_twoBody(s, mEta, mEta), phsp_twoBody(s, mEta, mEtap);
645 
646  Eigen::Array<fpcomplex, NCHANNELS, NCHANNELS> F = getPropagator(kMatrix, phaseSpace, adlerTerm);
647 
648  if(is_pole) { // pole
649  fpcomplex M = 0;
650  for(int i = 0; i < NCHANNELS; i++) {
651  fptype pole = couplings(i, pterm);
652  M += F(0, i) * pole;
653  }
654  return M / (POW2(pmasses(pterm)) - s);
655  } else { // prod
656  return F(0, pterm) * (1 - s0_prod) / (s - s0_prod);
657  }
658 }
659 #endif
660 
661 __device__ fptype phsp_FOCUS(fptype s, fptype m0, fptype m1) {
662  fptype mp = (m0 + m1);
663  fptype mm = (m0 - m1);
664  fptype a2 = (1.0 - mp * mp / s) * (1.0 - mm * mm / s);
665  return sqrt(a2);
666 }
667 
668 __device__ fpcomplex FOCUSFunction(fptype Mpair, fptype m1, fptype m2, unsigned int *indices) {
669  fptype s = POW2(Mpair);
670  unsigned int mod = GOOFIT_GET_INT(8);
671 
672  // mKPlus, mPiPlus, mEtap
673  constexpr fptype sNorm = mKPlus * mKPlus + mPiPlus * mPiPlus;
674  constexpr fptype s12 = 0.23;
675  constexpr fptype s32 = .27;
676 
677  fptype I12_adler = (s - s12) / sNorm;
678  fptype I32_adler = (s - s32) / sNorm;
679 
680  fptype rho1 = phsp_FOCUS(s, mKPlus, mPiPlus);
681  fptype rho2 = phsp_FOCUS(s, mKPlus, mEtap);
682 
683  fptype pmass = 1.7919;
684  Eigen::Array<fptype, 2, 1> coupling;
685  coupling << 0.31072, -0.02323;
686  fptype X = s / sNorm - 1;
687 
688  // constructKMatrix
689  Eigen::Array<fptype, 2, 2> kMatrix;
690 
691  kMatrix(0, 0) = coupling(0) * coupling(0) / (pmass - s) + 0.79299 - 0.15099 * X + 0.00811 * POW2(X);
692  kMatrix(1, 1) = coupling(1) * coupling(1) / (pmass - s) + 0.17054 - 0.0219 * X + 0.00085655 * POW2(X);
693  kMatrix(1, 0) = coupling(1) * coupling(0) / (pmass - s) + 0.15040 - 0.038266 * X + 0.0022596 * POW2(X);
694  kMatrix(0, 1) = coupling(0) * coupling(1) / (pmass - s) + 0.15040 - 0.038266 * X + 0.0022596 * POW2(X);
695 
696  fptype K11 = I12_adler * kMatrix(0, 0);
697  fptype K12 = I12_adler * kMatrix(0, 1);
698  fptype K22 = I12_adler * kMatrix(1, 1);
699 
700  fptype K32 = I32_adler * (-0.22147 + 0.026637 * X - 0.00092057 * POW2(X));
701 
702  fptype detK = K11 * K22 - K12 * K12;
703  fpcomplex del{1 - rho1 * rho2 * detK, -(rho1 * K11 + rho2 * K22)};
704 
705  fpcomplex T11{1., -rho2 * K22};
706  fpcomplex T22{1., -rho1 * K11};
707  fpcomplex T12{0., rho2 * K12};
708 
709  fpcomplex T32 = 1. / fpcomplex(1, -K32 * rho1);
710 
711  if(mod == static_cast<unsigned int>(Lineshapes::FOCUS::Mod::Kpi))
712  return fpcomplex(K11, -rho2 * detK) / del;
713  else if(mod == static_cast<unsigned int>(Lineshapes::FOCUS::Mod::KEta))
714  return K12 / del;
715  else /*if(mod==Lineshapes::FOCUS::Mod::I32)*/
716  return T32;
717 
718  return {0., 0.};
719 }
720 
721 __device__ resonance_function_ptr ptr_to_LS_ONE = LS_ONE;
722 __device__ resonance_function_ptr ptr_to_BW_DP4 = BW;
723 __device__ resonance_function_ptr ptr_to_lass = lass_MINT;
724 __device__ resonance_function_ptr ptr_to_glass3 = glass_MINT3;
725 __device__ resonance_function_ptr ptr_to_bugg_MINT = bugg_MINT;
726 __device__ resonance_function_ptr ptr_to_bugg_MINT3 = bugg_MINT3;
727 __device__ resonance_function_ptr ptr_to_SBW = SBW;
728 __device__ resonance_function_ptr ptr_to_NONRES_DP = nonres_DP;
729 __device__ resonance_function_ptr ptr_to_Flatte = Flatte_MINT;
730 __device__ resonance_function_ptr ptr_to_Spline = Spline_TDP;
731 __device__ resonance_function_ptr ptr_to_FOCUS = FOCUSFunction;
732 #if GOOFIT_KMATRIX
733 __device__ resonance_function_ptr ptr_to_kMatrix = kMatrixFunction;
734 #endif
735 
736 // This constructor is protected
737 Lineshape::Lineshape(std::string name)
738  : GooPdf(name) {
739  // Making room for index of decay-related constants. Assumption:
740  // These are mother mass and three daughter masses in that order.
741  // They will be registered by the object that uses this resonance,
742  // which will tell this object where to find them by calling setConstantIndex.
743 }
744 
745 std::vector<fptype> make_spline_curvatures(std::vector<Variable> vars, Lineshapes::spline_t SplineInfo) {
746  size_t size = std::get<2>(SplineInfo) - 2;
747  Eigen::Matrix<fptype, Eigen::Dynamic, Eigen::Dynamic> m(size, size);
748  for(size_t i = 0; i < size; i++) {
749  m(i, i) = 4;
750  if(i != size - 1) {
751  m(i, i + 1) = 1;
752  m(i + 1, 1) = 1;
753  }
754  }
755  m = m.inverse();
756 
757  Eigen::Matrix<fptype, Eigen::Dynamic, 1> L(size);
758  for(unsigned int i = 0; i < size; ++i)
759  L[i] = vars[i + 2].getValue() - 2 * vars[i + 1].getValue() + vars[i].getValue();
760 
761  auto mtv = m * L;
762  fptype spacing = (std::get<0>(SplineInfo) - std::get<1>(SplineInfo)) / std::get<2>(SplineInfo);
763 
764  std::vector<fptype> ret(vars.size(), 0);
765  for(unsigned int i = 0; i < size; ++i) {
766  ret.at(i + 1) = 6 * mtv(i) / POW2(spacing);
767  }
768  return ret;
769 }
770 
771 Lineshapes::GSpline::GSpline(std::string name,
772  Variable mass,
773  Variable width,
774  unsigned int L,
775  unsigned int Mpair,
776  FF FormFac,
777  fptype radius,
778  std::vector<Variable> AdditionalVars,
779  spline_t SplineInfo)
780  : Lineshape(name) {
781  GOOFIT_ADD_PARAM(2, mass, "mass");
782  GOOFIT_ADD_PARAM(3, width, "width");
783 
784  GOOFIT_ADD_INT(4, L, "L");
785  GOOFIT_ADD_INT(5, Mpair, "Mpair");
786 
787  GOOFIT_ADD_INT(6, enum_to_underlying(FormFac), "FormFac");
788 
789  GOOFIT_ADD_CONST(7, radius, "radius");
790 
791  if(std::get<2>(SplineInfo) != AdditionalVars.size())
792  throw GeneralError("bins {} != vars {}", std::get<2>(SplineInfo), AdditionalVars.size());
793  GOOFIT_ADD_CONST(8, std::get<0>(SplineInfo), "MinSpline");
794  GOOFIT_ADD_CONST(9, std::get<1>(SplineInfo), "MaxSpline");
795  GOOFIT_ADD_INT(10, std::get<2>(SplineInfo), "NSpline");
796 
797  int i = 11;
798  for(auto &par : AdditionalVars) {
799  GOOFIT_ADD_PARAM(i++, par, "Knot");
800  }
801 
802  std::vector<fptype> SplineCTerms = make_spline_curvatures(AdditionalVars, SplineInfo);
803 
804  for(auto &par : SplineCTerms) {
805  // Calc curve
806  GOOFIT_ADD_CONST(i++, par, "CKnot");
807  }
808 
809  GET_FUNCTION_ADDR(ptr_to_Spline);
810 
811  GOOFIT_FINALIZE_PDF;
812 }
813 
814 Lineshapes::GLASS::GLASS(std::string name,
815  Variable mass,
816  Variable width,
817  unsigned int L,
818  unsigned int Mpair,
819  FF FormFac,
820  fptype radius,
821  std::vector<Variable> AdditionalVars)
822  : Lineshape(name) {
823  if(5 != AdditionalVars.size()) {
824  throw GeneralError("It seems you forgot to provide the vector with the five necessary variables for GLASS, a, "
825  "r, phiF, phiR and F (in that order)");
826  }
827 
828  GOOFIT_ADD_PARAM(2, mass, "mass");
829  GOOFIT_ADD_PARAM(3, width, "width");
830 
831  GOOFIT_ADD_INT(4, L, "L");
832  GOOFIT_ADD_INT(5, Mpair, "Mpair");
833 
834  GOOFIT_ADD_INT(6, enum_to_underlying(FormFac), "FormFac");
835 
836  GOOFIT_ADD_CONST(7, radius, "radius");
837 
838  for(int i = 0; i < 5; i++) {
839  GOOFIT_ADD_PARAM(8 + i, AdditionalVars[i], "LassVars");
840  }
841 
842  GET_FUNCTION_ADDR(ptr_to_glass3);
843 
844  GOOFIT_FINALIZE_PDF;
845 }
846 
847 Lineshapes::One::One(
848  std::string name, Variable mass, Variable width, unsigned int L, unsigned int Mpair, FF FormFac, fptype radius)
849  : Lineshape(name) {
850  GOOFIT_ADD_PARAM(2, mass, "mass");
851  GOOFIT_ADD_PARAM(3, width, "width");
852 
853  GOOFIT_ADD_INT(4, L, "L");
854  GOOFIT_ADD_INT(5, Mpair, "Mpair");
855 
856  GOOFIT_ADD_INT(6, enum_to_underlying(FormFac), "FormFac");
857 
858  GOOFIT_ADD_CONST(7, radius, "radius");
859 
860  GET_FUNCTION_ADDR(ptr_to_LS_ONE);
861 
862  GOOFIT_FINALIZE_PDF;
863 }
864 
865 Lineshapes::RBW::RBW(
866  std::string name, Variable mass, Variable width, unsigned int L, unsigned int Mpair, FF FormFac, fptype radius)
867  : Lineshape(name) {
868  GOOFIT_ADD_PARAM(2, mass, "mass");
869  GOOFIT_ADD_PARAM(3, width, "width");
870 
871  GOOFIT_ADD_INT(4, L, "L");
872  GOOFIT_ADD_INT(5, Mpair, "Mpair");
873 
874  GOOFIT_ADD_INT(6, enum_to_underlying(FormFac), "FormFac");
875 
876  GOOFIT_ADD_CONST(7, radius, "radius");
877 
878  GET_FUNCTION_ADDR(ptr_to_BW_DP4);
879 
880  GOOFIT_FINALIZE_PDF;
881 }
882 
883 Lineshapes::LASS::LASS(
884  std::string name, Variable mass, Variable width, unsigned int L, unsigned int Mpair, FF FormFac, fptype radius)
885  : Lineshape(name) {
886  GOOFIT_ADD_PARAM(2, mass, "mass");
887  GOOFIT_ADD_PARAM(3, width, "width");
888 
889  GOOFIT_ADD_INT(4, L, "L");
890  GOOFIT_ADD_INT(5, Mpair, "Mpair");
891 
892  GOOFIT_ADD_INT(6, enum_to_underlying(FormFac), "FormFac");
893 
894  GOOFIT_ADD_CONST(7, radius, "radius");
895 
896  GET_FUNCTION_ADDR(ptr_to_lass);
897 
898  GOOFIT_FINALIZE_PDF;
899 }
900 
901 Lineshapes::NonRes::NonRes(
902  std::string name, Variable mass, Variable width, unsigned int L, unsigned int Mpair, FF FormFac, fptype radius)
903  : Lineshape(name) {
904  GOOFIT_ADD_PARAM(2, mass, "mass");
905  GOOFIT_ADD_PARAM(3, width, "width");
906 
907  GOOFIT_ADD_INT(4, L, "L");
908  GOOFIT_ADD_INT(5, Mpair, "Mpair");
909 
910  GOOFIT_ADD_INT(6, enum_to_underlying(FormFac), "FormFac");
911 
912  GOOFIT_ADD_CONST(7, radius, "radius");
913 
914  GET_FUNCTION_ADDR(ptr_to_NONRES_DP);
915 
916  GOOFIT_FINALIZE_PDF;
917 }
918 
919 Lineshapes::Bugg::Bugg(
920  std::string name, Variable mass, Variable width, unsigned int L, unsigned int Mpair, FF FormFac, fptype radius)
921  : Lineshape(name) {
922  GOOFIT_ADD_PARAM(2, mass, "mass");
923  GOOFIT_ADD_PARAM(3, width, "width");
924 
925  GOOFIT_ADD_INT(4, L, "L");
926  GOOFIT_ADD_INT(5, Mpair, "Mpair");
927 
928  GOOFIT_ADD_INT(6, enum_to_underlying(FormFac), "FormFac");
929 
930  GOOFIT_ADD_CONST(7, radius, "radius");
931 
932  GET_FUNCTION_ADDR(ptr_to_bugg_MINT);
933 
934  GOOFIT_FINALIZE_PDF;
935 }
936 
937 Lineshapes::Bugg3::Bugg3(
938  std::string name, Variable mass, Variable width, unsigned int L, unsigned int Mpair, FF FormFac, fptype radius)
939  : Lineshape(name) {
940  GOOFIT_ADD_PARAM(2, mass, "mass");
941  GOOFIT_ADD_PARAM(3, width, "width");
942 
943  GOOFIT_ADD_INT(4, L, "L");
944  GOOFIT_ADD_INT(5, Mpair, "Mpair");
945 
946  GOOFIT_ADD_INT(6, enum_to_underlying(FormFac), "FormFac");
947 
948  GOOFIT_ADD_CONST(7, radius, "radius");
949 
950  GET_FUNCTION_ADDR(ptr_to_bugg_MINT3);
951 
952  GOOFIT_FINALIZE_PDF;
953 }
954 
955 Lineshapes::Flatte::Flatte(
956  std::string name, Variable mass, Variable width, unsigned int L, unsigned int Mpair, FF FormFac, fptype radius)
957  : Lineshape(name) {
958  GOOFIT_ADD_PARAM(2, mass, "mass");
959  GOOFIT_ADD_PARAM(3, width, "width");
960 
961  GOOFIT_ADD_INT(4, L, "L");
962  GOOFIT_ADD_INT(5, Mpair, "Mpair");
963 
964  GOOFIT_ADD_INT(6, enum_to_underlying(FormFac), "FormFac");
965 
966  GOOFIT_ADD_CONST(7, radius, "radius");
967 
968  GET_FUNCTION_ADDR(ptr_to_Flatte);
969 
970  GOOFIT_FINALIZE_PDF;
971 }
972 
973 Lineshapes::SBW::SBW(
974  std::string name, Variable mass, Variable width, unsigned int L, unsigned int Mpair, FF FormFac, fptype radius)
975  : Lineshape(name) {
976  GOOFIT_ADD_PARAM(2, mass, "mass");
977  GOOFIT_ADD_PARAM(3, width, "width");
978 
979  GOOFIT_ADD_INT(4, L, "L");
980  GOOFIT_ADD_INT(5, Mpair, "Mpair");
981 
982  GOOFIT_ADD_INT(6, enum_to_underlying(FormFac), "FormFac");
983 
984  GOOFIT_ADD_CONST(7, radius, "radius");
985 
986  GET_FUNCTION_ADDR(ptr_to_SBW);
987 
988  GOOFIT_FINALIZE_PDF;
989 }
990 
991 Lineshapes::FOCUS::FOCUS(std::string name,
992  Mod mod,
993  Variable mass,
994  Variable width,
995  unsigned int L,
996  unsigned int Mpair,
997  FF FormFac,
998  fptype radius)
999  : Lineshape(name) {
1000  GOOFIT_ADD_PARAM(2, mass, "mass");
1001  GOOFIT_ADD_PARAM(3, width, "width");
1002 
1003  GOOFIT_ADD_INT(4, L, "L");
1004  GOOFIT_ADD_INT(5, Mpair, "Mpair");
1005 
1006  GOOFIT_ADD_INT(6, enum_to_underlying(FormFac), "FormFac");
1007 
1008  GOOFIT_ADD_CONST(7, radius, "radius");
1009 
1010  GOOFIT_ADD_INT(8, static_cast<unsigned int>(mod), "Lineshape modifier");
1011 
1012  GET_FUNCTION_ADDR(ptr_to_FOCUS);
1013 
1014  GOOFIT_FINALIZE_PDF;
1015 }
1016 
1017 #if GOOFIT_KMATRIX
1018 Lineshapes::kMatrix::kMatrix(std::string name,
1019  unsigned int pterm,
1020  bool is_pole,
1021  Variable sA0,
1022  Variable sA,
1023  Variable s0_prod,
1024  Variable s0_scatt,
1025  std::array<Variable, NCHANNELS> fscat,
1026  std::array<Variable, NPOLES *(NPOLES + 1)> poles,
1027  Variable mass,
1028  Variable width,
1029  unsigned int L,
1030  unsigned int Mpair,
1031  FF FormFac,
1032  fptype radius)
1033  : Lineshape(name) {
1034  GOOFIT_ADD_PARAM(2, mass, "mass");
1035  GOOFIT_ADD_PARAM(3, width, "width");
1036 
1037  GOOFIT_ADD_INT(4, L, "L");
1038  GOOFIT_ADD_INT(5, Mpair, "Mpair");
1039 
1040  GOOFIT_ADD_INT(6, enum_to_underlying(FormFac), "FormFac");
1041 
1042  GOOFIT_ADD_CONST(7, radius, "radius");
1043 
1044  GOOFIT_ADD_INT(8, pterm, "pTerm");
1045  GOOFIT_ADD_INT(9, is_pole ? 1 : 0, "Channel");
1046 
1047  GOOFIT_ADD_PARAM(10, sA0, "sA0");
1048  GOOFIT_ADD_PARAM(11, sA, "sA");
1049  GOOFIT_ADD_PARAM(12, s0_prod, "s0_prod");
1050  GOOFIT_ADD_PARAM(13, s0_scatt, "s0_scatt");
1051 
1052  for(int i = 0; i < NCHANNELS; i++) {
1053  GOOFIT_ADD_PARAM(14 + i, fscat.at(i), "fscat");
1054  }
1055 
1056  for(int i = 0; i < NPOLES * (NPOLES + 1); i++) {
1057  GOOFIT_ADD_PARAM(14 + NCHANNELS + i, poles.at(i), "Poles");
1058  }
1059 
1060  GET_FUNCTION_ADDR(ptr_to_kMatrix);
1061 
1062  GOOFIT_FINALIZE_PDF;
1063 }
1064 #endif
1065 
1066 Amplitude::Amplitude(std::string uniqueDecayStr,
1067  Variable ar,
1068  Variable ai,
1069  std::vector<Lineshape *> LS,
1070  std::vector<SpinFactor *> SF,
1071  unsigned int nPerm)
1072  : _uniqueDecayStr(std::move(uniqueDecayStr))
1073  , _ar(ar)
1074  , _ai(ai)
1075  , _SF(std::move(SF))
1076  , _LS(std::move(LS))
1077  , _nPerm(nPerm) {}
1078 
1079 bool Amplitude::operator==(const Amplitude &A) const {
1080  return _uniqueDecayStr == A._uniqueDecayStr && _ar == A._ar && _ai == A._ai && _LS == A._LS && _SF == A._SF
1081  and _nPerm == A._nPerm;
1082 }
1083 
1084 } // namespace GooFit
Generated on Sat Apr 21 2018 20:13:18 for GooFit by   doxygen 1.8.13