GooFit  v2.1.3
ThreeGaussResolution_Aux.cu
Go to the documentation of this file.
1 #include <cmath>
2 #include <goofit/PDFs/physics/ThreeGaussResolution_Aux.h>
3 
4 namespace GooFit {
5 
6 const fptype R1o6 = 1.0 / 6.0;
7 #define SQRTPIo2 (1.0 / M_2_SQRTPI)
8 #define SQRT1o2PI (sqrt(0.5 * M_1_PI))
9 
10 __device__ void gaussian(fptype &_P1,
11  fptype &_P2,
12  fptype &_P3,
13  fptype &_P4,
14  fptype _tau,
15  fptype adjTime,
16  fptype xmixing,
17  fptype ymixing,
18  fptype adjSigma) {
19  fptype _1oSqrtA = adjSigma * M_SQRT2;
20  fptype _1oSigma = 1 / adjSigma;
21  fptype _1o2SqrtA = 0.5 * _1oSqrtA;
22  fptype _1oSigma2 = _1oSigma * _1oSigma;
23  fptype _NormG = SQRT1o2PI * _1oSigma;
24 
25  fptype _C = 0.5 * adjTime * adjTime * _1oSigma2;
26  fptype _Bgn = -adjTime * _1oSigma2;
27 
28  fptype _Gamma = 1 / _tau;
29  fptype _B = _Gamma + _Bgn;
30 
31  fptype _u0 = _1o2SqrtA * _B;
32  fptype _u02 = _u0 * _u0;
33  fptype _F = _1oSqrtA * exp(-_C + _u02);
34 
35  fptype _Ig0 = SQRTPIo2 * erfc(_u0);
36  fptype _Ig1 = 0.5 * exp(-_u02);
37  fptype _Ig2 = _Ig1 * _u0 + 0.5 * _Ig0;
38  fptype _Ig3 = _Ig1 * (_u02 + 1);
39 
40  fptype _R = xmixing * _Gamma * _1oSqrtA;
41  fptype _R2 = _R * _R;
42  fptype _It0 = _F * _Ig0;
43  fptype _It1 = _F * _R * (_Ig1 - _u0 * _Ig0);
44  fptype _It2 = _F * _R2 * (_Ig2 - _u0 * _Ig1 * 2 + _u02 * _Ig0);
45  fptype _It3 = _F * _R2 * _R * (_Ig3 - _u0 * _Ig2 * 3 + _u02 * _Ig1 * 3 - _u0 * _u02 * _Ig0);
46 
47  // fptype _P0 = _NormG * _It0;
48  _P2 = _NormG * (_It0 - 0.5 * _It2);
49  _P4 = _NormG * (_It1 - R1o6 * _It3);
50 
51  fptype _u0py = _1o2SqrtA * (_B - ymixing * _Gamma);
52  fptype _u0my = _1o2SqrtA * (_B + ymixing * _Gamma);
53  fptype _Fpy = _1oSqrtA * exp(-_C + _u0py * _u0py);
54  fptype _Fmy = _1oSqrtA * exp(-_C + _u0my * _u0my);
55  fptype _Ipy = _Fpy * SQRTPIo2 * erfc(_u0py);
56  fptype _Imy = _Fmy * SQRTPIo2 * erfc(_u0my);
57  _P1 = _NormG * 0.5 * (_Ipy + _Imy);
58  _P3 = _NormG * 0.5 * (_Ipy - _Imy);
59 }
60 
61 __device__ fptype device_threegauss_resolution(fptype coshterm,
62  fptype costerm,
63  fptype sinhterm,
64  fptype sinterm,
65  fptype tau,
66  fptype dtime,
67  fptype xmixing,
68  fptype ymixing,
69  fptype sigma,
70  fptype *p,
71  unsigned int *indices) {
72  fptype coreFraction = RO_CACHE(p[RO_CACHE(indices[1])]);
73  // fptype tailFraction = p[indices[2]];
74  fptype tailFraction = (1 - coreFraction) * RO_CACHE(p[RO_CACHE(indices[2])]);
75  fptype outlFraction = 1 - coreFraction - tailFraction;
76  fptype coreBias = RO_CACHE(p[RO_CACHE(indices[3])]);
77  fptype coreScaleFactor = RO_CACHE(p[RO_CACHE(indices[4])]);
78  fptype tailBias = RO_CACHE(p[RO_CACHE(indices[5])]);
79  fptype tailScaleFactor = RO_CACHE(p[RO_CACHE(indices[6])]);
80  fptype outlBias = RO_CACHE(p[RO_CACHE(indices[7])]);
81  fptype outlScaleFactor = RO_CACHE(p[RO_CACHE(indices[8])]);
82 
83  fptype cp1 = 0;
84  fptype cp2 = 0;
85  fptype cp3 = 0;
86  fptype cp4 = 0;
87  gaussian(cp1, cp2, cp3, cp4, tau, dtime - coreBias * sigma, xmixing, ymixing, coreScaleFactor * sigma);
88  fptype tp1 = 0;
89  fptype tp2 = 0;
90  fptype tp3 = 0;
91  fptype tp4 = 0;
92  gaussian(tp1, tp2, tp3, tp4, tau, dtime - tailBias * sigma, xmixing, ymixing, tailScaleFactor * sigma);
93  fptype op1 = 0;
94  fptype op2 = 0;
95  fptype op3 = 0;
96  fptype op4 = 0;
97  gaussian(op1, op2, op3, op4, tau, dtime - outlBias * sigma, xmixing, ymixing, outlScaleFactor * sigma);
98 
99  fptype _P1 = coreFraction * cp1 + tailFraction * tp1 + outlFraction * op1;
100  fptype _P2 = coreFraction * cp2 + tailFraction * tp2 + outlFraction * op2;
101  fptype _P3 = coreFraction * cp3 + tailFraction * tp3 + outlFraction * op3;
102  fptype _P4 = coreFraction * cp4 + tailFraction * tp4 + outlFraction * op4;
103 
104  fptype ret = 0;
105  ret += coshterm * _P1;
106  ret += costerm * _P2;
107  ret -= 2 * sinhterm * _P3;
108  ret -= 2 * sinterm * _P4; // Notice sign difference wrt to Mikhail's code, because I have AB* and he has A*B.
109 
110  // cuPrintf("device_threegauss_resolution %f %f %f %f %f\n", coshterm, costerm, sinhterm, sinterm, dtime);
111  return ret;
112 }
113 
114 __device__ device_resfunction_ptr ptr_to_threegauss = device_threegauss_resolution;
115 
116 ThreeGaussResolution::ThreeGaussResolution(
117  Variable cf, Variable tf, Variable cb, Variable cs, Variable tb, Variable ts, Variable ob, Variable os)
118  : MixingTimeResolution()
119  , coreFraction(cf)
120  , tailFraction(tf)
121  , coreBias(cb)
122  , coreScaleFactor(cs)
123  , tailBias(tb)
124  , tailScaleFactor(ts)
125  , outBias(ob)
126  , outScaleFactor(os) {
127  GET_FUNCTION_ADDR(ptr_to_threegauss);
128  initIndex();
129 }
130 ThreeGaussResolution::~ThreeGaussResolution() = default;
131 
132 void ThreeGaussResolution::createParameters(std::vector<unsigned int> &pindices, PdfBase *dis) {
133  pindices.push_back(8);
134  pindices.push_back(dis->registerParameter(coreFraction));
135  pindices.push_back(dis->registerParameter(tailFraction));
136  pindices.push_back(dis->registerParameter(coreBias));
137  pindices.push_back(dis->registerParameter(coreScaleFactor));
138  pindices.push_back(dis->registerParameter(tailBias));
139  pindices.push_back(dis->registerParameter(tailScaleFactor));
140  pindices.push_back(dis->registerParameter(outBias));
141  pindices.push_back(dis->registerParameter(outScaleFactor));
142 }
143 
144 fptype ThreeGaussResolution::normalisation(
145  fptype di1, fptype di2, fptype di3, fptype di4, fptype tau, fptype xmixing, fptype ymixing) const {
146  // NB! In thesis notation, A_1 = (A + B), A_2 = (A - B).
147  // Here di1 = |A^2|, di2 = |B^2|, di3,4 = Re,Im(AB^*).
148  // Distinction between numerical subscribts and A,B is crucial
149  // for comparing thesis math to this math!
150 
151  fptype timeIntegralOne = tau / (1 - ymixing * ymixing);
152  fptype timeIntegralTwo = tau / (1 + xmixing * xmixing);
153  fptype timeIntegralThr = ymixing * timeIntegralOne;
154  fptype timeIntegralFou = xmixing * timeIntegralTwo;
155 
156  fptype ret = timeIntegralOne * (di1 + di2); // ~ |A|^2 + |B|^2
157  ret += timeIntegralTwo * (di1 - di2); // ~ Re(A_1 A_2^*)
158  ret -= 2 * timeIntegralThr * di3; // ~ |A|^2 - |B|^2
159  ret -= 2 * timeIntegralFou * di4; // ~ Im(A_1 A_2^*)
160 
161  return ret;
162 }
163 } // namespace GooFit
Generated on Sat Apr 21 2018 20:13:18 for GooFit by   doxygen 1.8.13